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ABSTRACT 


We  investigate  an  inverse  imaging  problem  for  TWI  (Through-the-wall  Imag¬ 
ing)  using  frequencies  between  500  GHz  and  1  THz.  Starting  from  first  principles,  this 
thesis  uses  Maxwells  equations  to  develop  a  model  for  the  transmission  Greens  func¬ 
tion.  This  simplified  model  is  then  used  in  a  Lippman-Schwinger  integral  equation  to 
predict  the  scattered  held  associated  with  interrogating  THz  waves.  We  investigate 
the  effects  of  wave  propagation  through  isotropic  media,  and  present  methods  for  cre¬ 
ating  images  from  the  scattered  held.  These  methods  are  examined  using  simulated 
data. 
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I.  INTRODUCTION 

Far  better  an  approximate  answer  to  the  right  question,  which  is  often 
vague,  than  an  exact  answer  to  the  wrong  question,  which  can  always  be  made 
precise.  John  W.  Tukey,  1962 

X-ray  computed  tomography  (CT)  was  created  in  the  late  1960’s  to  early 
1970’s.  Since  that  time,  this  imaging  system  has  come  to  be  recognized  as  one  of  the 
greatest  achievements  in  radiology  [Ref.  3].  From  this  point  in  time  to  the  present, 
solving  inverse  problems  and  ill-posed  problems  has  become  an  exciting  and  significant 
held  for  research.  As  in  reverse  engineering  approaches,  an  “inverse  problem”  is 
simply  one  which  attempts  to  attain  information  about  a  “source”  from  available 
information. 

Unlike  a  “direct  problem”  in  which  the  descriptive  parameters  of  the  problem 
are  known  in  advance  of  the  solution,  and  the  system  of  equations  can  be  solved  (sub¬ 
ject  to  initial  and  boundary  conditions),  an  inverse  problem  poses  a  special  challenge 
because  the  parameters  are  not  known  up-front. 

Suppose  we  want  to  measure  the  electromagnetic  held  at  a  given  position 
caused  by  a  signal  x ,  and  the  held  data  d  can  be  written  as 

d  =  {Ax}  (1.1) 

where  A  is  an  operator  that  describes  the  nature  of  the  measurement  system. 

Now  let  us  entertain  the  idea  that  the  object  x  E  Rn  is  mapped  into  the  data 
space  d  E  Rm  via  the  operator  A  which  would  be  a  mapping  of  the  object  from  the 
Hilbert  space  X  into  another  Hilbert  space  Y.  This  problem  is  considered  well-posed 
if  it  meets  the  conditions  in  Table  I: 

Conditions  1  and  2  tell  us  that  A  is  injective  and  spans  all  of  Y.  The  require¬ 
ment  of  condition  3  is  a  necessary  but  not  sufficient  condition  for  the  stability  of  the 
solution.  The  third  condition  is  motivated  by  the  fact  that  in  all  applications  the  data 
will  be  measured  quantities.  If  the  problem  is  well-posed,  the  error  propagation  is 
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1.  A  solution  exists 

2.  The  solution  is  unique 

3.  The  inverse  mapping  is  continuous 


Table  I.  Conditions  for  well-posed  problems 

controlled  by  the  condition  number.  For  example,  if  there  are  n  variables,  there  must 
be  n  independent  equations  [Ref.  30]  in  order  to  obtain  a  solution.  Therefore  the 
mapping  of  n  variables  by  the  system  of  equations  to  the  data  space  can  be  written 
as  Axn  —  d  where  A  is  a  matrix  defining  the  system  of  equations  and  xn  is  a  vector 
of  n  variables.  By  condition  1  and  2,  A”1  exists. 

Now  let’s  say  that  there  is  a  small  error  Ad  in  the  measured  data.  Then  the 
corresponding  uncertaintyAa:n  for  the  object  Axn  obeys 


I  Ax  , 


Xr 


<  cond(A) 


1M1 

lldll 


The  condition  number  is  defined  as 


cond(A)  =  1 1 A 1 1  || A 


-ii 


where  the  norm  1 1 A 1 1  is 


I  |Acc, 

=  sup  — ; — 

X„^0  II 

If  cond(A)  is  not  large,  then  the  problem  is  considered  to  be  well-conditioned, 
and  the  solution  will  be  stable  with  respect  to  small  variations  of  data.  If  the  condition 
number  is  large,  we  classify  the  problem  as  ill-conditioned.  The  difference  between  ill- 
posed  and  ill-conditioned  is  therefore  somewhat  vague.  Hadamard  defines  a  problem 
to  be  ill-posed  if  it  does  not  satisfy  all  three  conditions  of  Table  I.  An  ill-posed 
problem  is  one  where  an  inverse  does  not  exist  because  d  +  Ad  is  outside  the  range 
of  A.  Generally,  more  than  one  image  can  map  into  the  same  data,  or  small  changes 
in  the  data  can  result  in  large  changes  in  the  image. 
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Most  inverse  problems  are  ill-posed.  Therefore,  no  unique  answer  exists.  Let 
us  consider  the  Fredholm  integral  equation,  which  we  will  discuss  in  depth  later  in  this 
thesis.  In  real-world  problems,  the  integral  has  a  square-integrable  (Hilbert-Schmidt) 
kernel  and  is  of  the  form: 

k(x,  s)f(s)d.s  =  d(x)  a  <  x  <  b 

If  /  has  a  small  variation,  say  of  the  form  A /  =  esin(27uus),  uj  =  1, 2, 3, . . .,  e  G  R 
and  is  small.  Then  the  variation  of  data  will  be 

fb 

A d(x)  —  e  /  k(x,  s)  sm(2Trujs)ds 

J  a 

In  Appendix  A,  we  show  that  as  u  —>  oo  then  Ad  — ►  0.  Therefore,  the  ratio  ^ 
heads  to  oo  as  a;  heads  to  oo.  This  problem  is  ill-posed,  and  the  example  demonstrates 
that  Fredholm  integral  equations  can  be  sensitive  to  high-frequency  variations. 

Because  of  this  behavior  much  research  has  been  done  in  solving  these  types 
of  problems,  and  many  methods  and  techniques  have  been  developed.  This  thesis 
will  explore  some  of  these  techniques  and  create  a  simple  electromagnetic  scattering 
model  for  application  to  an  inverse  problem  for  high  frequency  imaging,  appropriate 
to  frequencies  ranging  from  10  GHz  -  1  THz. 

A.  BACKGROUND 

In  recent  years,  terahertz  technology  has  become  an  exciting  field.  Demands 
for  further  development  have  increased  many  fold  for  military,  security,  civilian  and 
medical  applications.  As  solid  state  physicists  are  coming  up  with  equipment  that  can 
provide  stable  terahertz  radiation,  and  as  detectors  are  getting  more  sophisticated 
in  measuring  these  data,  scientists  and  engineers  are  now  starting  to  focus  on  the 
terahertz  region  with  excitement.  This  is  largely  due  to  the  fact  that:  at  these 
frequencies,  there  is  great  material  penetration  capability;  unlike  X-rays,  the  radiation 
is  non-ionizing;  unlike  ultrasound,  you  can  image  without  contact;  and  penetration 


3 


depth  is  generally  greater  than  near-IR  radiation  [Ref.  8].  The  most  useful  region  of 
this  electromagnetic  spectrum  is  from  0.3  to  20  THzfRef.  16]. 

B.  SOME  PROPERTIES  IN  THE  GIGA-TERAHERTZ 
FREQUENCIES 

Frequencies  of  terahertz  (a  period  of  lps)  have  some  important  properties. 
Most  electrons  in  highly  excited  atomic  Rydberg  states  orbit  at  this  frequency.  Col¬ 
lisions  between  gas-phase  molecules  at  room  temperature  last  about  1  ps.  Super¬ 
conducting  energy  gaps  are  found  at  THz  frequencies.  Important  collective  modes 
of  proteins  are  known  to  vibrate  at  THz  frequencies  [Ref.  16].  These  are  just  a  few 
examples  of  terahertz  characteristics  in  the  physical  world.  This  region  lies  above 
the  frequency  range  for  most  electronics  but  below  the  range  of  optical  and  infrared 
equipment;  therefore,  it  has  been  difficult  to  explore  until  recent  developments  in 
solid  state  technology. 

However,  there  are  limitations  in  propagating  at  terahertz  frequencies  because 
they  are  highly  susceptible  to  atmospheric  absorption.  Only  a  small  window  within 
this  region  is  practicable. 

The  Terahertz  Imaging  Focal  Plane  Technology  (TIFT)  program  at 
DARPA  (Arlington,  VA)  defines  the  upper  limit  of  this  window  to  be  557 
GHz. 

There  are  also  certain  frequencies  beyond  this  optimum  band  that  can  provide  good 
propagation  behavior.  These  windows  are  centered  on  650GHz,  850GHz  and  at 
1.5THz.  Anything  above  1.5  THz  results  in  absorption  that  is  currently  too  high 
to  gather  any  meaningful  data[Ref.  8]. 

C.  MOTIVATION 

As  devices  provide  reliable  radiation  at  THz  frequencies,  there  is  also  a  need 
for  more  sophisticated  scattering  models.  There  are  many  models  already  in  place 
for  other  fields  such  as  radar,  and  we  have  had  much  success  in  extracting  superb 
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images  in  SAR,  X-ray  tomography,  sonograms,  etc.  If  these  powerful  methods  can  be 
carried  over  to  the  terahertz  range,  it  may  become  both  a  rewarding  and  a  worthwhile 
goal  for  development  in  the  area.  Recently,  through-the-wall  high  frequency  imaging 
has  become  a  promising  new  technology  that  addresses  many  civilian  and  military 
problems.  Much  of  the  practical  challenges  lie  in  the  reliability  of  the  system,  its 
portability,  and  its  ability  to  sustain  an  acceptable  level  of  performance. 

One  of  the  greatest  challenges  in  imaging  at  these  frequencies  is  that  the  system 
parameters  are  affected  by  the  structured  ambiguities  in  the  walls,  and  this  can  cause 
difficulty  in  measurement.  One  has  to  account  for  diffraction  and  refraction  effects, 
but  the  fact  that  we  can  image  with  relative  success  has  shown  that  this  area  has  a 
promising  future. 

There  are  other  factors  that  make  terahertz  radiation  an  attractive  area  for 
research.  Waves  at  this  frequency  can  penetrate  non-metallic  and  non-polar  media; 
moreover,  explosives,  chemical  agents,  and  many  different  biological  agents  return  a 
unique  spectra.  Furthermore,  this  radiation  is  safe  for  the  operator  of  the  device,  and 
is  harmless  to  the  biological  life  that  it  may  scan. 

Until  the  devices  that  generate  terahertz  signals  can  be  made  portable,  prac¬ 
tical  applications  are  still  over  the  horizon.  Academic  institutions  such  as  Rensse¬ 
laer  Polytechnic  Institute  have  demonstrated  the  reliabliity  of  these  sources  and  are 
currently  developing  a  more  portable  system  for  which  identification,  imaging,  and 
various  real  world  applications  can  be  used. 

However,  there  are  few  other  devices  that  are  already  in  place  that  are  small 
and  reliable.  Argonne  National  Laboratory  (ANL)  is  currently  investigating  whether 
such  penetration  capability  can  be  achieved  at  lower  frequencies  (in  the  10  —  500  Ghz 
region).  Such  waves  would  be  carried  through  a  wave  guide,  and  the  return  fields 
are  measured  via  a  waveguide.  ANL  demonstrated  that  these  waves  do  have  great 
penetrating  capabilities,  and  they  are  able  to  form  images.  But  imaging  performance 
at  these  frequencies  is  an  open  question. 
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Though  imaging  can  be  done,  most  methods  take  advantage  of  imaging  via 
bi-static  transmission  where  the  transmitter  is  on  one  end,  and  the  receiver  on  its 
opposite  end  so  that  the  image  can  be  done  via  intensity  strengths  of  field  going 
through  the  media  and  the  object  of  the  image.  However,  imaging  from  the  back- 
scattered  field  is  more  of  a  challenge  because  of  the  problems  encountered  not  only  in 
the  perturbation  caused  by  the  object  but  also  the  interactions  that  occurs  through 
the  transmission  medium. 

This  thesis  will  emphasize  the  modeling  requirements  using  Maxwell’s  equa¬ 
tions  and  will  develop  inversion  techniques  for  image  reconstruction.  We  will  first 
discuss  the  basic  Inverse  Problem.  Then  we  will  derive  (from  Maxwell’s  equations) 
the  theory  for  the  simple  model  used  in  this  thesis;  derive  the  Green’s  function  for 
each  region  in  space;  establish  the  scattered  field  equation  (perturbation  of  the  field 
due  to  the  object);  discuss  the  various  different  approximations  that  can  be  done  to 
establish  the  Fredholm  equations  for  the  scattered  field;  use  the  Born  approxima¬ 
tion  and  parametize  the  equations;  discuss  the  Singular  Value  Decomposition  and 
the  Least  Square  Method;  establish  the  problems  of  error  due  to  noise  and  discuss 
regularization;  and  finally  conduct  some  image  reconstruction  using  simulated  data. 
The  final  chapter  will  summarize  the  results  of  this  thesis  and  provide  some  ideas  for 
a  more  in-depth  model  for  image  reconstruction. 

D.  PRELIMINARY  DISCUSSION 
1.  Material  Effects 

For  our  model,  we  assume  that  the  wall  characteristics  are  known  with  a 
linear,  homogeneous,  and  isotropic  dielectric  constant.  There  are  certain  ranges  of 
Terahertz  that  are  not  very  useful.  For  example,  ranges  above  5  to  15  THz  frequencies 
have  highly  strong  phonon  absorption  in  virtually  all  crystalline  materials  [Ref.  11]. 
Innately  metallic  targets  cause  full  terahertz  reflection.  Skin,  due  to  the  high  water 
content,  absorbs  the  terahertz  rays,  and  most  of  the  energy  is  dissipated  in  the  first  100 
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microns  [Ref.  12].  One  must  balance  between  the  image  object  and  the  frequency  that 
would  provide  the  optimum  penetration  of  the  medium  while  providing  the  optimum 
reflection  of  the  object.  If  we  are  looking  for  metallic  material,  it  is  an  easier  case 
because  the  focus  would  then  be  penetrating  the  medium. 

2.  Atmospheric  Effects 

Another  important  consideration  is  how  feasible  is  it  to  transmit  the  electro¬ 
magnetic  waves  at  these  frequencies.  There  will  be  transmission  losses  in  the  field’s 
strength  through  the  material  and  reflection  of  the  target,  but  if  atmospheric  ab¬ 
sorption  is  too  strong,  one  may  just  as  well  abandon  the  idea  of  imaging  at  these 
frequencies.  Absorption  and  scattering  in  transmission  occur  largely  due  to  the  inter¬ 
actions  of  many  different  elements  in  the  atmosphere.  However,  absorption  is  mainly 
caused  by  oxygen  and  water  vapor:  the  greatest  amount  of  absorption  occurs  when 
the  frequencies  match  the  natural  resonances  of  oxygen  and  water  vapor.  As  the 
frequencies  deviate  away  from  the  natural  resonant  frequency  of  the  elements,  “win¬ 
dows”  (frequencies  at  which  absorption  is  minimal)  are  created.  These  are  the  areas 
in  which  we  are  especially  interested.  The  Liebe  model  gives  us  a  simple  view  of  the 
attenuation  and  dispersive  affects  in  the  atmosphere  [Ref.  14].  In  Figure  1,  we  can  see 
that  frequencies  between  0.7  THz  to  1  THz  have  some  great  areas  where  attenuation 
would  be  small  enough  that  we  can  neglect  absorption  effects. 

If  we  work  at  frequencies  <  1  terahertz  (such  as  0.8  THz),  we  have  almost 
no  dispersive  delay,  and  our  attenuation  is  also  quite  small  (<  O.lOdB/km).  Because 
atmospheric  affects  determine  the  type  of  devices  and  the  methods  of  gathering  data, 
we  have  to  first  determine  which  frequencies  would  give  us:  (1)  the  least  attenuation; 
(2)  the  least  dispersive  effects;  and  (3)  the  most  stable  and  common  atmospheric 
configuration.  In  Figure  2,  we  see  that  absorption  is  too  high  for  any  frequencies 
higher  than  1  THz,  and  it  will  remain  so  until  we  get  into  the  infrared  and  visible 
range  —  we  will  not  be  able  to  get  any  meaningful  data  in  this  region.  Therefore, 
most  of  the  research  and  work  should  be  within  this  limit. 
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Figure  1.  Attenuation  a  and  delay  (3  rates  for  three  fog  events  (w  = 
0.10,  0.25,  andQ.bg/ m 3)  added  to  a  saturated  air  (RH=100%)  sea  level  condition.  This 
Figure  is  taken  from  [Ref.  14] 
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Figure  2.  Atmospheric  Attenuation  from  [Ref.  15] 


E.  MODELING  CRITERIA 

Now  that  the  unexplored  electromagnetic  spectrum  is  starting  to  become  acces¬ 
sible,  the  areas  of  research  are  growing  without  bounds.  For  behind-the-wall  imaging, 
the  interest  for  us  is  to  develop  a  model  that  could  predict  held  information  using 
the  Green’s  function  which,  in  turn,  yields  the  system  (PSF:  Point  Spread  Function). 
The  utility  of  having  an  analytical  form  for  the  PSF  is  that  this  is  the  kernel  of  the 
Fredholm  Integral  equation.  Therefore,  knowing  the  nature  of  the  kernel  which  pro¬ 
vides  the  mapping  behavior  from  one  Hilbert  Space  to  the  other,  we  can  then  find 
the  source  (object)  of  the  perturbed  fields.  Since  we  can  measure  the  fields,  we  can 
solve  for  the  object  x  in  equation  1.1. 

To  solve  for  the  PSF,  we  will  make  the  following  assumptions: 


1.  The  waves  are  not  dispersive:  they  each  pass  through  the  medium  with  speed 

dw  _  „ 
dk  —  c 

2.  The  medium  is  linear  and  isotropic. 

3.  The  medium  is  non-magnetic. 

4.  The  medium  is  non-conductive. 
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Note  that  the  last  two  conditions  allow  us  to  simplify  Maxwell’s  equations. 
By  making  the  transmission  medium  non-conductive,  we  allow  the  index  of  refraction 
to  stay  real  and  no  damping  will  take  place.  This  is  a  rather  simple  model,  but  one 
that  can  be  adjusted  in  the  future  to  implement  this  very  real  effect.  Making  the 
medium  non-magnetic  will  make  the  slowing  of  the  wave  to  be  only  dependent  on  the 
dielectric  constant.  And  lastly  to  prevent  the  coupling  of  Maxwell’s  equations,  we  are 
going  to  make  the  medium  linear  and  isotropic. 
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II. 


THE  INVERSE  PROBLEM 


This  chapter  is  intended  to  be  a  basic  overview  of  the  inverse  problem  that 
this  thesis  will  cover.  We  will  review  the  basic  ideas  of  ill-posed  problems  and  regular¬ 
ization  methods  for  forming  stable  approximate  solutions.  We  will  confine  ourselves 
to  linear  equations  with  compact  operators  in  Hilbert  space.  We  will  base  our  recon¬ 
struction  algorithm  design  on  the  Singular  Value  Decomposition  and  the  Tikhonov 
regularization  method.  There  are  many  regularization  methods  for  which  the  inverse 
problem  for  this  imaging  criteria  could  be  done.  The  intention  of  this  thesis  would 
limit  the  scope  of  the  most  basic  Tikhonov  regularization  method  and  the  optimum 
methods  for  imaging  would  be  achieved  stochastically.  For  more  in-depth  discussion 
of  ill-posed  problems  and  different  regularization  techniques  to  enhance  imaging,  refer 
to  Tikhonov  and  Arsenin  [Ref.  26]. 


A.  CONCEPTS  OF  ILL-POSEDNESS  REVISITED 

In  the  last  chapter,  we  saw  that  a  well-posed  problem  must  meet  the  three 
conditions  listed  in  Table  I.  We  now  examine  the  data  model  equation  in  the  form  of 
a  Fredholm  integral  equation: 

A(xi,  s)f(s)ds  =  d(xi )  a  <  s  <  b  c  <  x\  <  X2  <  ■  ■  ■  <  xn  <  d 


The  interval  [c,d]  is  partitioned  into  N  parts  for  convenience  to  represent  discrete 
measurements.  This  integral  can  be  simplified  into  the  following  matrix  problem: 

(  m  \ 


d{xi)  =  <  lim  'V  A(xi,  a  +  j As)  f  (a  +  j As)  As 

I  M— »oo  ' 


where  As  =  ( b—a)/M  (II. 1) 


or 


where 


d  =  lim  Af 

M—>o o 
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d(x\) 

A(xi,a  +  As) 

A(xi,a  +  2As) 

A(xi,a  +  MA.s) 

d  = 

d(x2) 

,  A  = 

A{x  2,  a  +  As) 

A(x2,a  +  2As) 

A(x2,a  +  MAs) 

d{x  N) 

A(xn,  a  +  As) 

A(xn,  a  +  2As)  •• 

■  A(xn,  a  +  MAs) 

f(a  +  As) 
f(a  +  2  As) 

/(a  +  MAs) 


A  will  define  the  kernel  which  is  to  operate  on  /  E  X  where  X  belongs  to  the 
Hilbert  space.  From  the  three  conditions  of  Table  I,  we  say  that  the  mapping  of  A 
is  therefore  one-to-one  and  spans  all  of  the  vector  space  Y.  Because  we  are  dealing 
with  compact  operators  [Ref.  27],  Y  will  also  be  a  member  of  the  Hilbert  Space. 
This  tells  us  that  the  inverse  A~l  must  be  continuous  and  bounded.  A  well  behaved 
problem  requires  that  f  E  X  and  the  data  d  E  Y  are  well  controlled  by  the  norms 
of  the  vectors  in  X  and  Y.  But  in  most  real  life  applications,  we  do  not  have  this 
continuous  inverse  mapping,  which  suggests  the  contamination  of  the  data  space  with 
elements  not  belonging  to  A(X)  C  Y.  This  results  in  the  instability  of  the  solution 
for  f  which,  as  discussed  in  the  previous  chapter,  is  caused  by  the  condition  number. 
Colton  [Ref.  27]  proves  that  Af  ~  d  is  improperly  posed  if  /  E  X  is  not  of  finite 
dimension.  In  everything  we  encounter  in  the  real  world,  /  is  infinite  in  dimension, 
and  therefore,  ill-posedness  of  the  equation  II.  1  is  unavoidable. 

Straightforward  solutions  of  ill-posed  problems  often  result  in  non-physical 
answers.  For  this  reason,  one  must  be  very  careful  in  the  discretization  of  equation 
II.  1.  One  may  think  that  making  the  operator  A  more  refined  would  make  the 
solution  more  reliable,  but  the  opposite  is  true  [Ref.  27].  Therefore,  techniques  must 
be  introduced  to  construct  stable  solutions. 


As 
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B.  METHOD  OF  LEAST  SQUARES 

A  standard  technique  in  mathematical  and  statistical  modeling  is  to  find  a 
“least  squares”  curve  fit  to  a  set  of  data  points.  This  technique  was  first  developed 
independently  by  Carl  Friedrich  Gauss  and  A.  M.  Legendre.  Whenever  we  measure 
data,  we  are  going  to  get  errors  in  measurements  or  inaccuracies.  Therefore,  we 
cannot  expect  the  curve  to  perfectly  fit  through  all  the  data  points.  Instead  we  want 
the  best  approximation  in  the  sense  that  the  sum  of  squares  of  the  range  between 
the  data  point  and  its  corresponding  curve  is  minimized.  A  least  squares  problem  is 
generally  formulated  for  overdetermined  systems.  For  example,  suppose  x  E  and 
b  E  Rm  are  related  by  the  system  of  equations  b  —  Hx.  In  general,  we  may  or  may 
not  find  a  vector  x  E  R"  to  find  b  E  Mm.  We  can,  however,  find  a  vector  x'  whose 
transformation  by  H  would  result  in  the  minimum  distance  between  b  and  Hx! . 


H(X) 


Figure  3.  b  E  Rm  while  x'  E  X  and  X  = 
We  form  the  residual: 


r  —  b  —  Hx' 


and  we  minimize  the  distance 


||r||  =  ||6  —  Hx'\\ 

From  the  geometry,  we  see  that  the  minimum  distance  between  b  and  Hx'  has  to  be 
orthogonal  to  the  plane.  It  can  be  shown  by  using  Pythagorean  Law  that  this  will 
result  in  the  least  squares  solution  [Ref.  7].  From  Figure  3,  one  can  see  that  Hx'  is 
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the  projection  of  b  onto  the  H(X).  Due  to  orthogonality,  r  must  lie  in  the  nullspace 
of  H ,  where  a  null  space  is  defined  as  the  set  of  all  vectors  x  G  Rn  such  that  for  a 
given  operator  A ,  Ax  —  0.  In  set-builder’s  notation  we  have: 

Null(H)  =  {x  E  M  n  :  Hx  =  0}  (II.2) 

By  the  Fundamental  Subspaces  Theorem  [Ref.  29]  which  says  if  A  is  an  m  x  n 
matrix,  then  Null(A)  =  Range(AT)-L  and  Null(AT )  =  Range(A)1-  and  equation 
II. 2,  HTr  =  0.  Therefore,  to  solve  the  least  squares  problem  Hx  =  6,  we  must  solve 

r  =  b  —  Hx' 

Hlr  =  HT(b  -  Hx') 

0  =  H'b  -  H'Hx' 

H'b  =  H'Hx' 

x'  -  (HTH)~1HTb 

and  {HTH)  lHT  is  a  projection  matrix  of  b  (otherwise  known  as  the  Moore  Penrose 
Inverse).  We  will  apply  this  method  of  model  fitting  using  regularization  technique 
in  the  next  section. 

C.  BRIEF  DISCUSSION  OF  REGULARIZATION  METH¬ 
ODS 

In  our  previous  discussion,  we  saw  that  the  residual  falls  nicely  into  the 
nullspace  of  HT,  but  unfortunately  in  “real  life”  problems,  this  residual  will  not 
completely  fall  under  the  null  space  because  it  will  be  contaminated  by  some  other 
element  that  we  now  include  as  noise.  Methods  for  constructing  a  stable  approximate 
solution  of  an  ill-posed  problem  caused  by  the  effects  mentioned  are  called  regular¬ 
ization  methods.  If  error/noise  is  introduced  into  the  equation  d  =  Af  such  that 
d'  =  Af  +  n,  then  the  perturbation  would  be: 

\\d'  —  d\\  <  ||n||  (II. 3) 
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d  will  belong  in  the  range  of  A,  but  d!  may  or  rnaynot  be  in  the  range  of  A.  From 
functional  analysis,  we  introduce  a  linear  compact  operator  R  such  that: 

lim  RaAf  =  f  a  >  0  (II. 4) 

a— 

This  is  called  a  regularization  scheme  for  the  operator  A  [Ref.  27].  a  is  considered 
the  regularization  parameter.  For  example,  equation  II. 4  implies 

lim  Rad  — *•  A~ld 

a— >0 

As  a  goes  to  0,  Ra  converges  to  A-1. 

Colton  [Ref.  27]  points  out  that  the  regularization  schemes  for  compact  op¬ 
erators  are  not  uniformly  convergent.  Because  physical  objects  are  generally  infinite 
in  dimension,  RaA  cannot  converge  to  its  norm  as  a  — »  0.  For  this  reason,  we  must 
search  for  a  good  regularization  parameter.  The  regularization  scheme  approximates 
/  of  the  data  model  Af  =  d  with  a  regularized  object  f'  satisfying 

f  =  Rad'  (II.5) 


Including  measurements,  we  have: 

f  -  f  =  Rad'  ~  f 

f  ~  f  =  Rad!  -  Rad  +  Rad  -  f 
f  -  f  =  Rad'  -  Rad  +  RaAf  -  f 

f~f  =  Ra(d'  -d)  +  RaAf  -  f  (II. 6) 

where  d!  —  d  =  n,  the  (unknown)  measurement  noise.  By  the  triangle  inequality,  we 
get 

||  f  ~  /||  <  |H|  1 1  Ra  1 1  +  ||  RaAf  ~  f  ||  (II. 7) 

The  error  in  object  estimation  depends  on  two  terms  of  equation  II. 7  [Ref.  27]:  The 
first  term  demonstrates  how  the  error  is  propagated  by  the  measurement  limitation 
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and  the  noise  of  the  system;  the  second  term  shows  that  error  comes  from  the  inverse 
mapping  A-1  from  specific  elements  in  Y  mapped  by  A  to  X,  and  the  linear  operator 
R  that  maps  all  elements  in  Y  to  X.  We  desire  that  as  the  regularization  parameter 
a  goes  to  0,  so  will  the  term  \  \RaAf  —  f\\.  Unfortunately,  this  destabilizes  the  first 
term  of  eqn  II. 7.  Since  Ra  — >  A~  ,  ||i?a||  becomes  larger  (and  thus  the  condition 
number).  This  behavior  rapidly  increases  the  reconstruction  error.  The  first  term 
in  eqn  II. 7  therefore  requires  a  larger  regularization  parameter  for  stabilization,  and 
the  second  term  requires  the  parameter  to  go  to  0  for  accuracy.  Therefore,  a  method 
must  be  reached  to  optimize  the  right  hand  side  of  equation  II. 7. 

One  way  to  accomplish  this  is  to  make  a  choice  that  the  regularization  param¬ 
eter  a  be  dependent  on  the  noise  level,  that  is,  a  =  a(n ).  From  equation  II. 3,  we  see 
that  d!  — >  d  as  n  — >  0,  and  so 

Ra{n)d'  -r  A~ld  =  f 

There  are  many  approaches  to  accomplish  this  behavior.  In  practice,  a  has 
to  be  determined  via  a  priori  or  a  posteriori  criteria.  In  most  of  the  cases  for  TWI 
(Through-the- Wall-Imaging),  we  are  not  certain  what  we  are  looking  for,  which  makes 
a  priori  less  practical  than  the  a  posteriori  choice.  For  proof  of  concept,  this  the¬ 
sis  will  generate  data  with  point  scatterers  (series  of  Dirac  Delta  functions),  which 
we  would  then  use  the  analytic  form  of  the  kernel  to  determine  the  point  scatters 
themselves.  We  assume  we  have  some  prior  knowledge  of  the  model,  which  would 
not  likely  be  the  case  in  a  “real  world”  type  application.  We  will  not  cover  a  more 
detailed  Bayesian  statistical  analysis  of  the  model  in  this  thesis,  which  is  required  for 
a  posteriori  criteria,  but  we  leave  it  for  future  research  to  improve  on  the  model. 

D.  TIKHONOV  REGULARIZATION 

We  note  from  the  previous  section  that  as  the  regularization  parameter  goes 
to  0,  we  have  an  instability  of  1 1 Ra \  |  as  a  — »  0.  This  is  primarily  due  to  the  fact 
that  the  noise  contamination  of  the  data  is  related  to  the  (nonzero)  eigenvalues  of 
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Ra.  As  we  try  to  get  a  more  accurate  solution,  the  first  term  of  the  error  in  equation 
II. 4  grows  like  7^  (for  proof  refer  to  Colton  [Ref.  27]).  This  is  what  some  refer  to 
as  the  penalty  term  where  we  punish  the  equality  for  a  large  norm.  This  form  of 
regularization  is  due  to  Tikhonov  [Ref.  26]. 

There  are  a  few  things  we  can  do  to  minimize  this  effect.  We  seek  to  minimize 
equation  II. 7  in  the  least  square  sense. 

fa  =  argmin{||n||2  ||i?a||2  +  \\RaAf  -  /||2}  (II. 8) 

all/ 

Using  the  fact  that  a  is  a  function  of  noise,  we  can  introduce  this  regularization 
parameter  into  equation  II. 8.  Consequently  with  this  and  equations  II. 6  and  II. 7,  we 
can  arrive  at  the  Least  Square  solution  for  Ra  via  Newton’s  Method  [Ref.  27]. 

Ra  =  (al  +  AtA)~1At  (II. 9) 

Note  that  as  a  — >  0,  equation  II. 9  gives  us  the  Moore  Penrose  inverse,  the 
simple  least-squares  solution.  From  SVD,  we  know  that  A  can  be  written  in  terms  of 
right  singular  vectors  (V  matrix)  and  left  singular  vectors  ( U  matrix)  with  common 
singular  values  Uj  such  that 

A  =  USVT 

where  S'  is  a  diagonal  matrix  with  shared  singular  values  (see  Appendix  C).  Because 
U  and  V  matrices  are  formed  from  eigenvectors  of  symmetric  matrices  of  AAT  and 
AtA  respectively,  U  and  V  are  orthogonal  matrices. 

At  =  ( USVT)T  =  VSTUT 

Now  we  make  the  subsitution  of  A  and  AT  into  equation  II. 9: 

Ra  =  ( aVIVT  +  VSUTUSTVT)~ lVSTUT 

Because  U  and  V  matrices  are  formed  from  eigenvectors  of  symmetric  matrices  of 
AAt  and  ATA  respectively,  U  and  V  are  orthogonal  matrices.  Therefore,  we  get 
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the  following: 

Ra  =  ( aVIVT  +  VS2VT)~1VSTUT 
Ra  =  [V(al  +  S2)ET]  _1  VSTUT 
Ra  =  V(al  +  S'2)' lVTVSTUT 


E.  TRUNCATED  SVD  REGULARIZATION 

Another  way  to  control  the  ill-posed  behavior  is  to  simply  truncate  the  SVD 
representation:  pick  a  small  threshold  value  a  and  compare  this  value  with  a  given 
singular  value,  which  you  will  determine.  Once  this  threshold  is  reached  you  replace 
the  corresponding  Da  by  0.  This  method  is  called  the  truncation  filter.  There  are 
many  ways  to  determine  where  the  threshold  would  lie.  In  this  thesis,  we  will  look  for 
where  the  steepest  changes  in  singular  value  occurs  and  will  truncated  at  that  point, 
most  of  which  would  be  trial  and  error. 
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Both  methods — Truncation  filtering  and  Tikhonov  regularization — are  effec¬ 


tive  at  reducing  or  minimizing  the  effects  of  noise.  However,  neither  technique  can 
independently  offer  a  “best”  threshold  value  a. 
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III.  WAVES  IN  MEDIA 

A.  PRELIMINARY  DISCUSSION 

In  this  chapter,  we  will  set  the  foundations  for  applying  Maxwell’s  Equations 
to  predict  held  perturbations  in  space  located  behind  a  transmission  medium.  We 
will  work  primarily  with  the  time  independent  wave  equations,  particularly  for  the 
electric  held.  In  free-space,  the  usual  approach  applies  the  free  space  Green’s  function 
to  the  source  term  of  the  inhomogenous  electromagnetic  wave  equation.  This  source 
is  otherwise  known  as  the  scattering  potential. 

If  there  is  a  source  of  perturbation  that  occurs  in  region  III  in  Figure  4  and 
we  have  a  receiver  or  measuring  device  in  region  I,  generally  the  free-space  Green’s 
function  would  work  if  region  II  is  of  the  same  medium  of  region  I  and  III,  that  is 
free-space.  However,  if  region  II  is  a  wall,  the  free  space  Green’s  function  does  not 
satisfy  the  conditions  determined  by  the  boundaries  of  the  wall.  This  will  be  ex¬ 
plained  in  the  following  sections  where  we  break  the  spherical  waves  into  half  space 
expansions  and  apply  the  appropriate  boundary  conditions.  We  must  also  consider 
the  effects  of  refraction  and  diffraction  in  the  media.  The  Green’s  function  therefore 
needs  to  be  solved  either  analytically  or  be  measured.  In  real  world  situations,  this 
function  (also  known  as  the  System  Transfer  Function  (STF))  cannot  usually  be  de¬ 
termined  analytically  and  is  often  measured.  Such  measurements  are  limited  by  the 
accuracy  of  the  device.  As  discussed  in  chapter  2,  noise  contamination  can  result  in 
catastrophic  image  reconstruction  artifacts  when  null-space  noise  components  con¬ 
taminate  the  data.  Furthermore,  the  reconstruction  becomes  more  unstable  as  the 
condition  number  increases.  Robust  and  powerful  techniques  of  regularization  will  be 
generally  required. 

In  the  next  chapter,  we  will  create  a  simple  model  where  we  approximate  this 
transfer  function  by  imposing  the  appropriate  boundary  conditions.  These  results  will 
then  be  applied  (in  the  spatial  frequency  regime)  to  the  development  of  an  imaging 
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algorithm. 


B.  PARAMETERS 

Electromagnetic  waves  through  matter  are  governed  by  three  properties  of 
the  material:  the  permittivity  e,  the  permeability  p  and  the  conductivity  a.  In  our 
model,  we  assume  the  material  to  be  non-magnetic  with  constant  permittivity  and 
no  conductivity,  which  means  the  index  of  refraction  is  going  to  be  constant  and 
real-valued.  We  will  discuss  some  of  the  effects  of  the  model  and  shortcomings  due 

to  this  simplification. 

a.  Maxwell’s  Equations 

Maxwell’s  equations  for  stationary  media  are  of  the  following  form  [Ref. 

6]: 


V  •  E 


V  •  B 


Vx£ 

VxB 


Ptotal  Ptrue  VP 

eo  eo 

0 

dB 

~~di 

( .  dP  „  dE\ 

p0  \  3true  +  +  V  X  M  +  60  —  I 


(III.l) 

(HI-2) 

(HI-3) 

(HI-4) 


where  ptotal  is  the  total  charge  density  composed  of  the  charges  that  are  free  to  move 
Ptme  and  those  charges  that  are  bounded  Pbounded  —  —  V  ■  P,  and  P  is  the  polarization 
of  the  medium.  jtrue  is  considered  the  true  or  the  free  current.  And  if  we  introduce 


I 


II 


III 


Figure  4.  Region 
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the  field  vectors  D  and  H  by 


D  =  e0E  +  P 

H  =  —-M 
P  o 


then  the  Maxwell  equations  simplify  to 


^  '  J-)  Ptrue 


V  •  B 


V  x  E 

V  x  H 


0 

dB 

"1 n 

dD 

3  true  3  oT 


(HI-5) 
(III. 6) 
(HI-7) 
(III.8) 


1.  Anisotropic  Media 

For  anisotropic  media,  the  permittivity  and/or  permeability  is  not  a  scalar- 
valued  parameter.  Since  most  materials  in  the  walls  are  non-magnetic,  we  can  safely 
conclude  that  permeability  will  remain  isotropic.  However,  most  “real”  problems 
may  have  anisotropic  permittivity.  Therefore  D  cannot  be  represented  by  a  simple 
constant  of  proportionality  to  E.  Instead  a  tensor  equation  is  necessary  so  that  [Ref. 
32] 

D  =  tensor(e)  ■  E  (III. 9) 

where 

E3 x  ^XxEx  T  C xyEy  -\-  CXzE z 

E)y  CyXE x  T  ^yyEy  "b  ^yzE z 

E) z  ^zxEx  T"  tzyEy  T  CZzE z 


and 


tens  or  (e) 


^ xx  £ xy  £ xz 

£yx  ^yy  ^ yz 

£ zx  £ zy  ^-zz 
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Since  D  may  not  be  parallel  to  E,  the  permittivity  tensor  is  a  function  of  nine 
numbers.  If  we  are  considering  a  homogeneous  media,  then  e.y  will  be  independent 
of  position.  Otherwise,  they  too  can  be  varying.  For  most  of  the  common  cases,  we 
expect  the  permittivity  tensor  to  be  symmetric  [Ref.  32].  This  means  that  only  6 
different  constants  are  needed  to  define  the  3x3  matrix.  Furthermore,  because  the 
matrix  is  symmetric,  its  eigenvectors  will  be  orthogonal.  Therefore,  if  you  choose  the 
coordinate  systems  such  that  it  is  aligned  to  these  eigenvectors,  you  can  reduce  the 
tensor  to  a  diagonal  tensor  whose  elements  are  its  eigenvalues  such  that 


tensor'  (e) 


4  0  0 

0  e'y  0 

0  0 


(III.10) 


This  coordinate  system  is  called  the  principal  system.  When  the  electric  field  lies 
along  one  of  these  principal  axis  (eigenvector),  then  D  ||  E,  and  thus  the  mathematical 
analysis  can  be  simplified. 

For  simplicity,  our  model  will  neglect  this  effect  (that  exists  in  materials  such 
as  crystals).  We  will  restrict  ourselves  to  a  linear  media  where  the  polarization  of  the 
media  is  directly  proportional  to  the  field. 


D  =  e0E  +  eo\E  =  e0(l  +  x)E  =  e0erE  =  eE  (III. 11) 

where  x  is  the  electric  susceptibility,  eo  is  the  permittivity  of  free  space,  er  is  the 
relative  permittivity,  and  e  is  the  dielectric  constant. 

2.  Conductive  Materials 

In  a  conductive  media,  which  largely  obeys  Ohm's  Law,  and  assuming  the 
material  is  linear,  we  write 

j  =  aE  =  —D  (III. 12) 

e 

By  the  equation  of  continuity 

V-i  +  |  =  0  (III. 13) 
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we  see  that 


=  0 

3  stationary 

SO  that  V  •  j stationary  =  0 

Substituting  equation  III.  12  into  this  last  result  yields 

•  edj  =  • 

J  '  ^  i  3  stationary 

3  3  stationary  T  j 0®  E  (III.  14) 

where  Trei  =  -  is  the  characteristic  or  the  relaxation  time  [Ref.  6]. 

Taking  the  curl  of  equation  III. 7,  we  obtain  the  following 

d 

Vx(Vx£)  =  --(VxB)  (III. 15) 

We  also  restrict  our  model  to  be  non-magnetic.  Most  wall  materials  are  indeed  non¬ 
magnetic,  and  we  are  safe  to  employ  this  approximation.  Therefore,  the  relative 
permeability  fir  ~  1  and  so  the  material  permeability  is  approximately  the  magnetic 
constant  in  vacuum  //  ~  yu0-  This  allows  us  to  write  B  =  ji^H.  With  Equation  III. 8, 
Equation  III.  15  then  can  be  simplified  to 


V-tj 


3  + 


dD\ 
~dt  ) 

dD 


dt 


a  /  f)F\ 

V (V  •  E)  -  V2E  =  (jtrue  +  e—J  (III. 16) 

For  most  cases,  we  will  be  dealing  with  charge  free  regions.  Hence  V  •  E  =  0. 
Now  for  conductive  material,  we  may  substitute  Equation  III. 12  into  III. 16  for  further 
simplification.  We  arrive  at  the  following  with  the  given  assumptions: 

d2  E  dE 

W'2E  ~  =  0  (IIL17) 
Harmonic  solutions  are  of  the  form: 
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E(x,  y,  z,  t )  =  E(x,  y,  z)ue  lU}t 


From  this  point,  we  shall  work  with  time  independent  wave  equations,  and 
will  refer  all  vectors  in  the  frequency  domain,  and  Equation  III. 17  becomes: 


V  E  y^eu/2  E  -(-  iy0aujE  —  0 
V2E  +  f1  +  !M)^2E  =  0 


(111. 18) 

(111. 19) 


Since  e  =  eoer,  c  = 


and  n 


■sJEr,  we  can  simplify  the  above  equation  to 


v2e+  ( 

2  ia  \ 

7l2  2 

is — ) 

—ujE 

\ 

c2 

V2£  + 

(i  +  A 

V  eo; . 

)  n2k2E 

(III. 20) 


where  oj/k  =  c.  With  the  relaxation  time  derived  earlier  rrei  =  e/a,  we  get 


V2E  +  1  + 


Trel^ 


n2k2E  =  0 


(III. 21) 


From  Equation  III.  18  we  define  a  complex  wave  number  k  =  k  +  in  so  that 


k2  —  yo  eu2  +  iy^auj 


(III. 22) 


and 


k 

K 


1/2 


+  1 


Note  that  for  all  pure  metals,  the  relaxation  time  is  ~  10-14  [Ref.  6],  which  is  on  the 
same  order  as  the  collision  time  of  the  electrons.  At  the  frequency  that  we  propose 
to  use  (500GHz  to  ITHz),  the  relaxation  time  will  be  of  the  order  jpJ  «  102  if  the 
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material  is  fully  conductive.  This  means  that  the  diffusion  part  of  the  equation  would 
dominate  and  therefore  the  attenuation  would  be  quite  high.  The  penetration  depth 
governed  by  d  =  1/k  will  be  quite  short,  and  the  waves  will  not  be  able  to  penetrate. 
Usually  for  a  good  dielectric  the  relaxation  time  is  of  the  order  of  days  [Ref.  33], 
and  hence  the  imaginary  term  of  Equation  III. 21  vanishes,  and  thus  there  will  be  no 
damping  in  the  media.  For  such  frequencies,  we  maybe  able  to  generalize  to  certain 
conductive  materials  with  larger  relaxation  times,  but  such  an  analysis  would  require 
a  large  database  of  material  effects  and  is  beyond  the  scope  of  the  thesis.  For  this 
reason,  the  model  also  assumes  that  the  walls  are  non-conductive  and  we  will  deal 
with  just  the  wave  equation  portion  of  Equation  III. 17 

3.  Frequency  Dependence  and  the  Dispersion  Relation 

Dispersive  media  have  values  of  yu,  e,  or  a  that  are  functions  of  frequency.  It 
can  be  shown  that  for  a  dielectric  e  as  a  function  of  frequency  is  well-modeled  by[Ref. 
33]: 


e(w) 

e\uj) 

e» 


e'(co)  —  ie"(uj) 

tQuliul-uj2) 

60  (u>2  -  u2)2  +  a2uj2 

e0  uj^uja 

(lo2  —  luq)2  +  oi2ijj2 


where  a  =  1/r  and  r  is  the  collision  times  of  the  electrons,  is  the  resonant  fre¬ 


quency. 


u2  =  -  plasma  frequency, 

N  are  the  number  of  dipoles  per  unit  volume,  and  m  is  the  mass  of  the  electron.  The 
real  part  e'(o;)  defines  the  refractive  index  n(ui )  =  yV (w)/e o,  and  the  imaginary  part 
e"(cu)  is  the  loss  tangent  of  the  material  tan#(w)  =  e"/e'(w),  and  it  is  related  to  the 
attenuation  constant  (absorption  coefficient)  of  the  electromagnetic  wave  [Ref.  33]. 
In  Figure  5,  we  see  that  around  the  resonant  frequency  luo,  the  dielectric  behaves  in 
a  strange  manner  and  the  material  absorption  becomes  quite  high.  Therefore,  one 
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has  to  be  careful  that  the  radiation  frequency  we  transmit  would  not  be  in  resonance 
with  the  dielectric,  and  further  material  analysis  might  be  required. 


Figure  5.  Real  and  Imaginary  parts  of  dielectric  constant.  Figure  taken  from  [Ref. 
33] 

For  metals,  charges  are  unbounded,  so  ujq  =  0,  and  from  Ohm’s  Law,  we  find 
[Ref.  33] 

a(u)  =  e°Up  (III. 23) 

In  “real”  phenomena,  we  know  that  from  optics  that  n  is  a  function  of  frequency 

(See  Figure  6).  Therefore,  it  can  be  seen  that  different  wavelengths  of  light  refract 

differently.  This  is  largely  due  to  dispersion  dtu/dk  /  constant.  Because  waves  of 

different  frequency  travel  at  different  speeds  in  a  dispersive  medium,  one  would  have 

to  find  methods  and  ways  to  isolate  these  frequencies  for  application.  The  waves  as 

a  whole  or  the  “envelope”  of  the  waves  travel  at  group  velocity: 

d  oj 
"9  =  dk 

Because  the  velocity  depends  on  the  frequency  and  oj  =  u)(k),  the  field  is  of 
the  form: 

/OO  /»oo 

A(k)eiikx~umdk  =  /  A(k)eit(kx/t_w(k))dk  (III.24) 

-oo  J  —  OO 

The  amplitude  A(k)  is  related  to  the  initial  condition  of  ip(x,  0)  by: 

/OO 

A(k)eikx  dk  (III. 25) 

-OO 
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Analysis  of  the  system  typically  requires  a  numerical  approach  (such  as  the  the 
Method  of  Stationary  Phase  [Ref.  5]).  The  motivation  for  this  is  that  we  know  that 
for  all  wave  numbers  k,  wave  component  disperse  for  any  initial  conditions.  All  the 
waves  will  travel  with  the  group  velocity,  and  all  of  the  energy  maybe  concentrated 
in  one  wave  number  [Ref.  5].  For  this  reason,  using  the  Method  of  Stationary  phase 
can  provide  us  a  very  good  prediction  of  the  system.  These  are  all  physical  realities 
we  will  have  to  face  in  dispersive  mediums,  for  which  case  a  large  database  of  various 
different  material  will  be  required.  It  is  beyond  the  scope  of  this  thesis  to  analyze  the 
effects  of  the  waves  for  each  different  media.  For  this  reason,  this  thesis  will  set  the 
criteria  that  we  are  working  with  non-dispersive  mediums. 

C.  FOUNDATIONS  OF  THE  SCALAR  THEORY 

From  our  Maxwell  equations  listed  in  Equations  (III. 5),  (III. 6),  (III. 7),  (III. 8), 
we  get  our  wave  equations  by  taking  the  curl  of  Equation  III. 7  and  Equation  III. 8. 
Since  the  vector  wave  equation  is  obeyed  by  both  E  and  H ,  an  identical  scalar  wave 
equation  is  obeyed  by  all  components  of  those  vectors  provided  we  are  in  a  cartesian 


Figure  6.  Index  of  Refraction  vs  Wavelength.  This  Figure  is  taken  from  [Ref.  31] 
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coordinate  system.  Thus,  the  component  Ex  will,  for  example,  obey  the  equation: 


V2EX  +  7i2k2Ex  =  0 

Therefore,  it  is  possible  to  summarize  the  behavior  of  all  components  of  E  and  H 
through  scalar  equations  of  the  form: 

V2^  +  n2k2ip  =  0  (III. 26) 

If  we  take  from  Equations  (III. 7)  and  (III. 8)  and  compare  component  by  com¬ 
ponent,  we  will  get  6  equations,  and  if  we  orient  the  coordinate  system  so  that  one 
component  vanishes,  then  the  equations  will  simplify  to  two  equations  with  the  E 
and  B  held.  These  equations  can  be  uncoupled  and  reduced  to  a  single  scalar  wave 
equation,  which  describes  either  a  Transverse  Magnetic  Polarization  (TM)  or  Trans¬ 
verse  Electric  Polarization  (TE).  When,  all  the  medium  parameters  e,  a,  and  yu  are 
independent  of  one  of  the  space  dimensions,  then  the  corresponding  electromagnetic 
held  components  are  also  independent  of  that  same  dimension.  Then  one  can  write 
the  vector  Maxwell  equations  as  a  superposition  of  the  TE  and  TM  helds. 

If  there  is  variation  of  permitivity  as  in  the  case  of  anisotropic  media,  the  wave 
equation  becomes 


V(E  ■  V(ln  e))  +  V2E  +  nV#  =  0  (III.27) 

In  this  case,  the  extra  term  couples  the  x,  y  and  z  components.  The  scalar  represen¬ 
tation  would  fail  in  the  conductor  case  where  the  boundary  conditions  are  that  there 
is  a  jump  of  a  in  the  electric  held  perpendicular  to  the  surface,  and  a  jump  with  the 
magnitude  of  the  free  current  in  the  magnetic  held  parallel  to  the  surface.  In  these 
cases,  the  components  of  the  Electric  and  Magnetic  helds  will  be  further  restricted, 
and  would  require  a  more  robust  model  such  as  the  Stratton-Chu  model  for  vector 
wave  equations  [Ref.  27]. 
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If  the  waves  are  spherical  in  nature,  one  can  still  formulate  the  problem  as 
a  scalar  wave  in  terms  of  the  Debeye  potentials  [Ref.  6].  In  either  case,  our  model 
criteria  gives  us  the  flexibility  for  using  the  scalar  model.  Although  there  are  faults  in 
using  the  scalar  model,  our  model  criteria  allow  us  to  decouple  the  Maxwell  equations 
and  reasonably  capture  the  essential  behavior  of  the  wave  propagation. 

Furthermore,  we  have  to  impose  the  conditions  at  infinity.  In  the  vector  case, 
the  electric  and  magnetic  fields  must  go  to  zero  at  infinity,  and  thus  must  satisfy  the 
Silver- Miiller  conditons: 


lim  (H  x  x  —  rE)  =  0  (III. 28) 

r— xx) 

lim  (E  x  x  +  rH)  =  0  (III. 29) 

r— »oo 

where  r  =  \x\  and  where  the  limit  is  assumed  to  hold  uniformly  in  all  directions. 
The  corresponding  radiation  condition  for  Cartesian  components  is  the  Sommerfeld 
radiation  condition 

f)  F1 

lim  r(— - ikE)  =  0  (III. 30) 

r— xx)  Of 

For  more  complete  discussion  and  proof,  refer  to  [Ref.  27]. 

For  the  reasons  mentioned  above,  we  shall  attack  the  problem  for  the  scalar 
case  with  the  Sommerfeld  conditions  imposed  at  infinity. 
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IV.  THE  INHOMOGENEOUS  WAVE 
EQUATION/TRANSFER  FUNCTION 


To  develop  the  transfer  function,  consider  a  field  that  emanates  towards  an 
infinite  flat  slab  (Figure  7).  We  take  this  wave  to  be  directed  in  the  positive  2 
direction,  and  y  axis  is  positive  downwards,  and  x  is  positive  outwards.  The  isotropic 
media  has  an  index  of  refraction  nm  that  is  independent  of  frequency,  and  the  location 
of  the  source  is  at  r0  =  —  z^k. 

For  a  held  with  source  s(r,t),  the  wave  equation  becomes 

where  £  is  in  the  time  domain  for  E.  The  source  function  can  be  written  in  terms  of 
the  Fourier  integral: 


Figure  7.  Point  Source  emanating  waves  through  a  constant  dielectric  medium  (In¬ 
finitely  large  slab) 
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s(r,t)=  /  S(f,r)ei2^  df 


Thus  the  time  independent  wave  equations  reduces  to 


V2E  +  n2k2E  =  S(f,r)  (IV.2) 

We  can  construct  the  solution  of  Equation  IV.2  by  superposition  of  unit  point  solu¬ 
tions  corresponding  to  a  source  at  the  point  r o  =  —zok,  and  we  placed  this  source  on 
the  z  axis  for  reasons  to  be  discussed  below. 

The  Green’s  function  due  to  a  source  at  vq  will  satisfy  equation: 


(V2  +  k2n2(z))G(x,  r0)  =  8(x  +  r0) 
where  k  =  2tt / X,  X  is  the  wavelength  in  free-space 

5(x  +  r0)  =  S(x)8(y)S(z  +  z0) 


and 


1  for 
n(z)  =  {  n  for 
1  for 

Therefore,  the  particular  solution  for  E  is 


z  <  0 
0  <  z  <  L 
L  <  z 


(IV.3) 


(IV.4) 


(IV.5) 


E=  [  G(x,  r0)S(f,  r)d3r0 
J  Vo 

and,  in  the  time  domain,  we  have 


£{r,t) 


E(r,  f)e~i2nft df 


(IV.6) 


(IV.7) 


Since  the  media  is  non-conductive,  the  wave  number  remains  real  valued. 
Therefore,  there  will  be  no  dampening  term,  and  the  media  will  become  a  wave  guide 
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composed  of  TE  and  TM  waves.  Using  Bracewell’s  notation,  we  take  the  Fourier 
transform  of  Equation  IV.  3  in  the  x  and  y  coordinate  to  obtain 

/OO  /»oo 

/  G(x,  y,  z,  r0)e-‘(2^+2^dxdy  (IV.8) 

•oo  J  —  OO 

/oo  /»oo 

/  G(/x,  fy ,  A  r0)e*(2^+2^dfxdfy  (IV.9) 

•oo  J  — oo 

This  transform  will  reduce  the  three  dimensional  Laplacian  to  a  single  ODE.  We 
are  allowed  to  make  this  substitution  because  there  is  no  slowing  or  changing  of  the 
media  in  the  x  and  y  direction:  the  media  remains  constant  in  these  directions,  and 
the  index  of  refraction  only  changes  in  the  z  direction.  Application  of  Equation  IV. 3 
yields 


p\  2 

[-(2irfx)2-(2TTfy)2+—  ]G(fxJy,z,ro)  +  k2n2{z)G(fxJy,z,r0)  =  5{z+z0)  (IV.10) 


Define 


F  =  2nfxi  +  2  irfyj 

so  that  (FI2  =  (27 r/x)2  +  (2tt fy)2. 

The  wave  number  is  related  as: 


Ax  Ay  Az 


=  (27 Tfx)2  +  (2vr  fy)2  +  (27 rfzf 


(IV.ll) 


(IV.12) 


Consequently,  the  Bracewell  notation  chooses  the  spatial  frequency  to  obey  fx  = 
1  /K  ,  fy  =  1/Ay  ,  fz  =  1/AZ  and  we  write  G  =  G(fx,fy,z,  r0) 

When  2  /  —  zo,  Equation  IV. 10  is  the  homogeneous  wave  equation  with  solu¬ 
tion 


G 


iy/k?n?{z)-\F\*z  +  b 

V  V  V 

- V -  - 

right  going  wave 


■iyj  k2n2  (z)—\F \2  z 

v  y 

left  going  wave 


(IV.  13) 
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This  solution  must  satisfy  the  boundary  conditions,  which  will  determine  the  coef¬ 
ficients  of  each  component  wave  as  they  propagate  through  the  medium.  The  wave 
and  its  derivatives  must  be  continuous  at  any  boundaries.  Consequently,  the  following 
conditions  must  be  true  (see  Figure  8): 


Y 

b- 

bm 

a 

(0,0, -Z_0) 

1 

am 

aB 

t 

A  “ 

n  m 

A  B 

G 

G 

G 

L 

Figure  8.  Cross  sectional  view  of  the  wavefront 


- 1  Am  | 

\z= 0  —  ^  \z= 0 

dG~ 

dGm  i 

dz 

8* 

II 

O 

qd  ! 

|  Z=L  =  GB\z—l 

dGm 

dz 

ddB 

z=L  dz 

(IV.  14) 


(IV.15) 


Gb 


_  c 

Z——ZO  ^ 


dG~ 


\Z=  —  ZQ 


dz 


dGb 


Z=  —  ZQ 


dz 


=  l 


(IV.  16) 


Z  =  -ZQ 


Imposing  the  Sommerfeld  conditions  and  using  the  fact  that  the  wave  is  outgo¬ 


ing  from  the  source,  the  coefficient  ab  of  the  right  going  component  (for  z  <  — £q),  and 
the  coefficient  bB  of  the  left  going  component  for  (z  >  L )  must  vanish  (see  Figure  8). 
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The  waves  are  traveling  to  the  left  when  z  <  —zq  and  are  traveling  to  the  right  when 
z  >  L.  The  boundary  conditions  are  such  that  the  Green’s  function  is  continuous 
with  continuous  derivatives  —  except  at  the  source.  Since  the  waves  are  continuous, 
the  Green’s  function  has  to  be  continuous.  But  its  first  derivative  is  not  very  obvious. 
Let  us  integrate  Equation  IV.  10  both  sides. 

G(fx,fy,z,r0)dz+ 1 (k2n2(z))  — (2vrfx)2  — (27rfy)2)G(fx,fy,z,r0)dz  =  J  S( z+z0)dz 

(IV.  17) 

The  right  hand  side  will  become  the  Heaviside  Function.  The  left  hand  side  would  be 
composed  of  the  first  derivative  of  the  Green’s  function  and  an  integral  of  the  Green’s 
function. 

f)  r 

■gzG(fx,  fy,  z,  r0)  +  J  (k2n2(z)  -  ( 2nfx )2  -  (2tt fy)2)G(fx,  fy,  z,  rQ) dz  =  H(z  +  z0) 

(IV.18) 

Since,  the  Green’s  function  is  continuous,  the  middle  term  will  vanish  at  any 
given  boundary.  However,  the  first  derivatives  will  become  continuous  as  long  the 
boundaries  are  z  <  —Zo  or  z  >  —zo-  As  the  first  derivative  crosses  the  source  (where 
the  Heaviside  Function  transitions  from  0  — >  1,  the  difference  across  will  have  a  jump 
of  1).  And  imposing  these  boundary  conditions  specified  above,  we  arrive  to  the 
following: 

Let 


At  z  =  0 


ut  =  n(z  <  0) 
nm  =  n( 0  <  z  <  L) 
ub  =  n(L  <  z) 


a~  +b~  =am  +  bm 


(IV.19) 
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^jk2n\  -  \F\2a~  -  <Jk2n2T  -  \F\2  = 

am^k2n2m  -\F\2  -  bm^/k2n2m  -  \F\2 


1  1 

a 

1  1 

am 

y/Nr  —  \/  Nt 

b~ 

VNm  —y/Nm 

bm 

Ml  M2 


where  Nt  =  k2nB  —  |_Fj2  and  Nm  =  A’2?z^  —  l-Fl2- 
At  z  =  L,  we  require 

gmeiy/k2n^n-\F\2L  e~iy/ k2nll-\F\'2 L  _  a_B gyj 'k2n2B-\F\2L 


(IV.20) 


(IV.21) 


(IV.22) 


am \fk2n2m  -  \F\2eiyJ^nl,-\F?L _ h™ ^/k2n2m  -  \F\2e~Wk2n2m-\F\2L 
=  aB ij k2n2B  -  \F\2ei^nB~\F\2L 


(IV.23) 


gi\/  Nml '  g  i'/  NmL 

am 

^^eiVN^L  Ze~iVW™L 

bm 

M3 


(IV.24) 


gi y/NsL  g— iV  iVsL 

"  " 

y/]y^'g*VA eA  g —iy/NsL 

0 

M  4 


where  Nm  =  A’2n^  —  \F\2  and  NB  =  k2n2B  —  (-F)2 
And  at  z  =  —zq  ,  we  require 


^>eiyJk2n\-\F\2ZQ  _  a~  e~iyJ k2ri^-\F\2z0  _|_  eiy/k2n^-\F\2  z0 


(IV.25) 


Now  we  apply  the  jump  condition  (see  Equation  (IV.  16))  to  obtain 
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(IV.26) 


a 


pi  1 2g ~iyj k2n'^,  —  \F\2zo 
-b~  ^Jk2n2T  -  |F|2eiVfc2«T-l^l2^ 

=  -bbJk2n2T-  \p\2eiVk^T-\^\2zo  +  i 


e~iy/NTz0  p^/Nrzo 

a 

y/N^e~i'jl^Z0  --v/iVye*^20 

b~ 

M  5 


(IV.27) 


giV  Ntzq  Q—iVNrzo 

—  sj  NTei'^'Z0  e~i'^'z° 

v- - ... - 

M  6 

where  iVy  =  k2n §,  —  |  ^  | 2 .  We  note  that  all  the  matrices  are  invertible,  and  therefore, 
we  can  have  a  closed  form  solution  for  the  coefficients  of  the  Green’s  function. 

To  solve  for  these  coefficients,  we  first  simplify  our  results  to  read 


Mx 

a 

=  m2 

am 

b~ 

bm 

M3 

am 

=  Mi 

"  aB  " 

bm 

0 

a 

bb 

Mh 

=  m6 

+ 

b~ 

0 

(IV.28) 


(IV.29) 

(IV.30) 


where  the  reflection  and  transmission  amplitudes  can  be  expressed  in  terms  of  a  , 
the  amplitude  of  the  incident  wave: 


„B 

a 

=  m1-1m2m3-1M4 

a 

b~ 

' - ^ - ' 

B 

0 

(IV.  31) 
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or 


B 


a  =  Bija 
b  -  =  -®21  -- 


B 


b  =  B2-\a 

R  d 

aB  = 


B 


1,1 


B 


1.1 


After  much  algebra,  we  obtain  the  reflectance: 


(IV.32) 

(IV.33) 


r=b^=  UNb  +  -v/A^X-y/Ar  -  vgQ  +  {-y/7h  +  ^)(y^Vr  +  VN^)ei2^L 
a~  {VNb  +  VAQ(-v/Ar  +  V^NQ  +  (-y^  +  v^Xv^  - 

(IV.34) 

The  Reflection  coefficient  then  is  R  =  r  *r. 

Similarly,  the  transmittance  is 


as  _  4y/N^^/NZe-iVW*L 

er  ~  (\/At  +  v^)(x/iVB  +  VN^)e~^L  +  (v^Vt  -  y/N^)e^L 

(IV.35) 


The  Transmission  is  T  =  t*t. 

We  can  also  solve  for  the  other  coefficients.  Inside  the  medium,  we  have 


am 

=  M^Ma 

'  aB  " 

bm 

c 

0 

(IV.36) 


am  =  ChlaB  =  Cl>1  um  „b  _  CVi 


Si,i 


=  C2)1a"  = 


B 


(IV.37) 


1,1 


a 


m 


a 


_ 2y^Vr(y^W:  +  Vgyjj) _ 

+  VKn){VN^  +  \/iVm)  +  (- -  VKn)e'^L 

(IV.38) 


6m 

a- 


_ 2-v/Ar(-y/A^  -  vgg) _ 

(v^Vb  +  VKnXVNr  +  V^m)  +  (- VNt  +  V^mXV^k  -  V/A 

(IV.39) 
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Finally,  at  the  source,  we  obtain 


bb 

0 


Me1 


^M5  Mf 1 M2  M3 1 M4 


0 


Writing  X  =  M6  1M5M1  1M2-M3  and 


(IV.40) 


yields 


S' 


0 

1 

ei'/Nr zo  (1+^/jVr) 

1 

1 

(IV.41) 


b6  -  XhlaB  -  Shl  =  (Xhl/Bhl)a-  -  Sltl  (IV.42) 

and 

0  =  X2^aB  -  S2,i  (IV.43) 

We  can  simplify  these  results  somewhat.  Outside  the  medium,  we  assume 
y/Nr  =  VX~b  =  VN.  We  are  interested  in  the  amplitudes  of  the  transmitted  wave 
(i aB ),  the  incident  wave  (a-)  and  the  reflected  wave  ( b~ ). 

In  the  next  chapter,  we  will  show  that  we  can  consider  the  Green’s  function 
to  be  shift-invariant  so  that  G(x,  r0)  =  G[x  —  Xq ,y  —  yo ,  z,  Zo)  =  G(r  —  r0). 

e  iVNz0 

a  =^vW 


(. N  -  Nm)  +  {Nm  -  N)ei2Vw^L  eiVNzo 

(VKn  +  s/N )2  -  (y/N^  -  \fN)2ei2'/w™L  2 y/N 


(IV.45) 


aB  =  t*  a  = 


4y/N^N^e~iVWL 


iy/~N  zq 


(IV.46) 


(y/N  +  -  (y/N/,  -  /N)2ei'^L  2y/N 

We  note  that  for  z  <  —z0,  is  free  space.  Since  the  free  space  Green’s  function  is 
a  solution  to  the  Helmholtz  equation  and  satisfies  the  radiation  condition  at  infinity, 
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it  must  be  the  solution  in  this  region.  We  will  show  this  for  —Zq  <  z  <  0  in  the  next 
chapter.  For  the  regions  of  interest,  we  have  the  following: 

Zo  <  z  <  0  right-going 
left-going 

(IV.47) 


G(fx,fy,z,r0 )  = 


eiVNz0  i  y/k2-\F\2z 
2  VN  C 

(■ N-Nm)+(Nrn-N)ei2'/N™L  e»yCVzQ  ; 

(\/N^+y/N)2-(VN^-y/N)'2ei2'/N^L  2 Vn  ^ 


42-|F|2 


G(fx,  fy,  z,  t*0)  =  a'VV'"2*2-^!2*  right  going  wave  (IV.51) 

G(fx,fy,z,r0)  =  bme-^n2k2~\F\^  left  going  wave  (IV.52) 
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V. 


THE  ANGULAR  SPECTRUM 


A.  WEYL’S  INTEGRAL  AND  SPHERICAL  WAVES 

In  this  chapter,  we  will  show  that  the  right-going  wave  of  Equation  IV. 47  is 
the  planar  expansion  of  the  spherical  wave  by  introducing  the  Weyl’s  integral.  We 
will  also  discuss  the  properties  of  the  Inverse  Transform  of  the  Green’s  function  in  the 
spatial  domain,  and  demonstrate  the  shift-invariant  form  of  the  Green’s  Function. 

What  makes  the  problem  very  difficult  is  the  reflection  and  refraction  of  a 
spherical  wave  on  a  planar  media  because  the  wave  has  spherical  symmetry  and  the 
boundary  is  a  plane.  For  this  reason,  we  want  a  way  to  expand  the  spherical  waves 
in  terms  of  planes.  In  equation  IV.47,  we  will  first  take  the  case  for  the  right-going 
wave.  The  Inverse  Fourier  Transform  defined  as: 


G(x,  y,  z,  r0)=  /  G(fx,  fy,  z,  r0)e^*+^dfxdfy 

J  — oo  J  — oo 

H.  Weyl’s  expansion  of  the  spherical  waves  into  planar  waves  is  [Ref.  24]: 


jk\r-r  |  jj,.  roo  poo  i 


It*  —  r'  |  2tt 


'  — OO  J  — OO 


_ ^Mx-x^+syiy-y^+s^z-z']]^ 

s7  v 


where  r  =  xi  +  yj  +  zk,  r’  =  x'i  +  y'j  +  z'k  and 


(V.l) 


(V.2) 


sz  =  +Jl-  s2x-  s2y  when  +  s2y  <  1 


sz  =  +iJ s2  +  s2  —  1  when  s2x  +  s2  >  1 
Therefore,  the  Inverse  Fourier  Transform  for  the  right  going  wave  becomes 


G(x,y,z,r0)  = 


poo  ei^k2-(2irfx)2-(27rfy)2(z+z0) 

l-oo  2ky/l  -  (27T fx/k)2  -  (2irfy/k)2 


ei(2nfxX+2nfyy)d^d^  (yg) 


Under  the  change  of  variables: 
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_  271 -fx 


dsr  — 


„  _  27 rjy 

* y  k 


27rdfx  j  _  2?rdfy 

k  Uby —  k 


Equation  V.3  becomes: 


/•oo  noo  ei^/k2-k2 s2-k2 s2(z+z0) ei(ksxx+ksyy)  /  fc  \  /  fc 

G{sx,sy,z,r0)=  I  I  - ______ - /_W_]dsxdsy 


’  —  oo  J  — oo 


2VX  ~  sl~  sl 

k  roo  no o  l-s2-s2{z+z0)+ksxx+ksyy) 


8i r2 


'  — oo  »/  — OO 


V1  -  sx  -  s 


dsxdsy 


and,  substituting  equation  V.4  into  equation  V.2,  we  arrive  at  the  following: 


G(x,y,z,r0)  =  - 


1  eik\r-r0\ 


■IAtt  \r  —  r0\ 

Operating  at  far  field,  this  Green’s  function  would  then  be  simplified  to 


eik\r-r0\  eikr 


0-ikr-T  o 


?’47r|r  —  r0 1  iAirr 

where  f  is  the  unit  vector  pointing  in  the  r  direction  [Ref.  25]. 


(V.4) 


(V.5) 


(V.6) 


B.  GREEN’S  FUNCTION  IN  MEDIA 

We  make  a  coordinate  transformation  of  Equation  IV. 47  and  Equation  IV. 48 
for  reasons  to  be  explained  in  the  following  sections.  Use  the  relationship  in  equation 
IV.  12,  and  substitute  the  spatial  frequencies  in  for  the  wave  number  such  that  we  get 
the  following  (see  Figure  9): 


2tt  fx  =  k  sin  6  cos  (j)  2nfy  =  k  sin  6  sin  <j)  (V.7) 

27 rfz  =  k  cos  6 

In  equation  IV. 47,  we  notice  that  e e lVk2-\F\2z  -s  t].ie  p]anar  expansion 
of  the  spherical  wave.  We  use  this  observation  to  see  some  physical  similarities  of 
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Figure  9.  Coordinate  Axis 


the  other  waves.  The  reflection  of  the  wave  in  Equation  IV. 47  can  be  written  as 
e~l^k2^^'\2z .  And  so  the  returning  wave  is  also  spherical  in  nature: 


R  = 


(■ N  -  Nm)  +  (Nm  -  N)ei2^L 
+  y/N )2  -  ( y/N. ~  -  y/N)^^L 


(V.  8) 


y/N^  =  \jn2mk2  -  (2 irfx)2  -  (2vr fy)2  y/N  =  \jk2  -  (2 irfx)2  -  (2vr fy)2  (V.9) 

From  equation  V.7,  we  rewrite  equation  V.9  as 

\f  Nm  —  kyj n ^  —  sin2  6  VN  —  k  \J\  —  sin2  6  —  k  cos  9  (V.10) 

Substituting  equation  V.10  into  V.8,  we  get  the  following: 

(1  -  n2)  +  (n2  -  1  )ei2fcWn2— sin20 

ItyU  j  —  . 

n2  —  sin2  0  +  -\/l  —  sin2  6)2  —  (yj n2  —  sin2  9  —  \/l  —  sin2  g^ei2kLVn2~sin2 e 

(V.ll) 

(1  -  n2)  +  (n2  -  i ynkLyJrV-Avte 
( \] n2  —  sin2  9  +  cos  9)2  —  (\/n2  —  sin2  6  —  cos  0)2ei2fcL^/n2-sin2e 
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The  Inverse  Fourier  Transform  defined  in  equation  V.l,  with  the  following 


change  of  variables: 


fx 


—  sin  6  cos  (j) 

Z7T 


fy 


—  sin  9  sin  6 
2tt 


so  that 


dfx 

dfx 

dd 

d<j> 

dfy 

dfy 

d6 

drj> 

cos  9  cos  (j)  —  ^  sin  9  sin  < 
cos  9  sin  (j)  ^  sin  9  cos  f 


The  spatial  reflected  Green’s  function  becomes 


—  1  cos  9  sin  9 

2i r 


G(x,y,z,r0)=  1 1  G(9,  r0)e<(te™®  “Basins  Bin*)  (  )  cos  9  sin  9d9d</>  (V.12) 

'  '  '27 r  1 


where  G(6,z,r0)  =  R(6)  e^c°7  e  lkz^ 

It  is  easy  to  see  that,  under  the  coordinate  transformation,  Equation  IV. 47 
can  be  made  into  Equation  V.12,  but  the  limits  of  integration  are  not  very  obvious. 
Because  we  are  summing  the  contributions  of  all  the  waves  from  various  propagation 
directions,  we  can  see  that  </>  varies  from  0  to  2tt.  The  projection  of  the  propagation 
vector  onto  the  x  —  y  plane  is  symmetric  for  any  given  angle  9  normal  to  the  me¬ 
dia.  However,  the  9  coordinate  requires  more  careful  analysis.  In  its  cartesian  form 
(Equation  V.l),  the  Inverse  Fourier  transform  suggests  that  there  will  be  an  area  of 
integration  where  fx  and  fy,  the  projected  components  of  the  propagation  vector  onto 
the  x  —  y  plane,  will  be  large  enough  that  the  argument  of  the  exponential  of  z  in  the 
frequency  domain  of  the  Green’s  function  is  no  longer  complex-valued  but  instead, 
becomes  exponentially  decaying  in  the  z  direction.  In  the  following  section,  we  shall 
discuss  this  effect,  and  discuss  how  9  will  have  to  be  adjusted  to  include  this  effect. 


C.  ANGLE  OF  REFLECTION/INCIDENCE 

Because  the  Green’s  functions  are  continuous  at  the  boundaries,  the  total  field 
at  the  interface  between  any  two  media  must  obey  incident  wave  +  reflected  wave  = 
transmitted  wave  for  all  the  regions  in  space  where  there  is  a  change  in  the  index  of 
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refraction.  Since  this  requirement  must  hold  true  for  all  time,  it  must  be  frequency 
independent.  Therefore 


lo(z  <  0)  =  <u(0  <  z  <  L)  —  uj(L  <  z )  (V.13) 

Now,  for  the  reflection  at  a  surface,  this  requirement  becomes 


c(z  <  0) bright— going  c{z  0)fcjey^_ 


going 


k 


right— going 


=  k, 


left— going 


( V.  14) 


Because  the  waves  must  be  continuous  at  the  interface,  the  projection  of  the  prop¬ 
agation  vector  onto  the  x  —  y  plane  will  be  the  same  at  the  boundaries (z  —  0  and 
z  =  L).  Therefore 


kx  =  2irfx,  ky  =  2Trfy,  and  k sin  9 

must  the  same  at  z  =  0  and  z  =  L  for  all  the  waves  (incident,  reflected,  transmitted). 
From  Equation  IV. 47,  a  single  wave  propagating  to  the  right  from  the  source  can  be 
written  as 

Gincident  =  G(fx,  fy,  z,  royW*x+^fyy)  (V.15) 


A\rN  zq  , - =— 

_  ® _ giy  k2-\F\2z ei(2nfxx+2nfyy) 

2  y/N 

and  the  reflected  wave  (the  left  going  wave)  is  written  as 


(V.16) 


G 


reflected 


(N  -  1 Vm)  +  (Nm  -  N)ea^-L  e^° 


(VKn  +  Xn)2  -  (s/Kn  -  ■jNYei^N~L  2 XN 


(V.17) 


Comparing  the  exponents,  we  see  that  only  the  ^-component  of  the  propagation  vector 


is  different  by  a  negative  sign  with  the  same  magnitude.  And  seeing  that  the  k  must 
be  equal,  we  conclude 

^k2  -  \F\2  =  yjk2  -  |F2|  (V.18) 

" - - - '  0. - v - ' 

incident  reflected 
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k  cos  6inc  =  k  cos  0ref  (V.  19) 

Oinc  =  Oref  (V.20) 

This  is  the  familiar  Law  of  Reflection. 

D.  SNELL’S  LAW 

For  transmission,  we  can  verify  Snell’s  Law  in  the  same  manner:  The  frequen¬ 
cies  and  the  projection  of  the  propagation  vector  at  the  x  —  y  plane  need  to  be  the 
same  due  to  the  boundary  conditions  discussed  earlier.  This  requirement  means: 


u>(z  <  0)  =  tu(0  <  z  <  L) 


ckrnr  —  A 


n 


trans 


nkjrtr  =  k 


trans 


If  we  measure  9trans  from  the  normal  (see  Figure  10),  we  can  see  from  the 
geometry  that 


sin  0 trans  ^ inc  Sin  0 incj Tlhi ■ 

Tt  sin  0t rans  sin  Sine 


(V.21) 


This  is  Snell’s  Law.  Equation  V.21  suggests  that  there  is  a  limit  to  the  angle  for 
transmission  because  the  sine  function  has  a  range  between  —1  and  1.  But  if  we  need 
to  integrate  further  into  infinity,  we  cannot  be  restricted  to  real  angles.  Stratton  [Ref. 
34]  introduces  the  complex  angle  of  incidence  through  which  the  sine  functions  will 
become  hyperbolic  functions.  This  approach  can  be  used  to  account  for  the  evanescing 
waves,  and  is  implicit  in  the  Inverse  Fourier  Transform  of  Equations  (IV. 47),  (IV.48), 
(IV.51),(IV.52).  In  the  Cartesian  form,  we  can  see  that  fz  — >  too  as  fx  —>  Too  or 
fy  -»•  ±oo. 

To  see  this  in  the  spherical  transformation,  we  examine  at  the  ^-component  of 
Equation  V.7:  cos  0  =  =  -\A  (27rR)  (2nfy)  .  j£  combinat;ion  0f  ari(]  js 

larger  than  k2,  then 
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Figure  10.  Snell’s  Law 


cos  9 


W^fxf  +  (27 Tfy)2 

k 


k 2 


roo 


(V.22) 


as  either  fx  or  fy  approaches  to  infinity.  From  the  geometry  in  Figure  9,  9  must 
go  from  0  to  7r/2,  but  which  branch  of  9  should  be  used  in  the  complex  plane?  A 
complex  angle  of  incidence  indicates  that  6  should  be  of  the  form  6  —  a  +  i/3,  and  so 
by  trig  identity: 


cos($)  =  cos(a  +  i/3 )  =  cos(a)  cos  (i/3)  —  sin(o:)  sin  (i/3)  (V.23) 

=  cos(a)  cosh(/3)  —  tsinh(/3)  sin(a)  (V.24) 

Our  allowable  angles  6  <  tt/2  and,  in  order  for  Equation  V.22  to  be  true,  f3 
must  go  towards  —  oo  since  sinh(/3)  is  an  odd  function.  Therefore,  we  choose  a  path 
for  which  9  =  0  — >  |  —  too. 


E.  EVANESCING  WAVES 

We  consider  the  region  —  zq  <  z  <  0.  From  Equations  (IV.47)  and  (IV. 51),  we 
consider  the  reflection  and  the  transmission  for  a  single  wave  front. 


G(x,y,z,r0)  =  A(fxJy)e~i^k2~{2nfx)2~{27Tfy)2{z+zo)ei{2nUx+2nfyy) 
G(x,y,z,r0 )  =  ^'(A,  A)eiV/^M20V MHf(^+^)ei(Hx+Hs) 

where  A(fx,fy)  and  A'(fx,  fy)  can  be  complex  or  real- valued  amplitudes. 


(V.25) 

(V.26) 
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Consider  the  reflected  wave:  the  propagation  vector  of  the  reflected  wave  must 
be  parallel  to  the  x—y  plane  before  the  wave  can  start  to  evanesce.  Once  this  condition 
exists,  we  can  see  that  the  exponent  in  the  z  direction  becomes  real,  and  the  only 
imaginary  portion  is  the  argument  2tt fxx  +  2Trfyy  =  F  ■  xi  +  yj.  Therefore,  the  wave 
propagation  is  parallel  to  the  media  before  it  can  start  to  evanesce.  But  from  a  larger 
index  of  refraction  to  air  (which  is  the  case  at  0  <  z  <  L)  where  the  rays  bend  away 
from  the  normal,  one  can  see  from  Equations  (IV. 48)  and  (IV.51)  that  there  is  an 
angle  of  incidence  will  result  in  a  transmission  angle  that  is  parallel  to  the  surface. 
This  is  the  condition  for  total  internal  reflection,  and  the  angle  is  the  critical  angle. 
From  this  view  point  as  the  z  component  of  the  incident  wave  inside  the  medium 
decreases,  the  angle  of  incidence  increases,  and  thus  the  transmitted  wave  will  start 
to  evanesce.  As  either  fx  or  fy  increases  further,  the  conditions  are  met  where  the 
reflected  waves  will  also  evanesce  in  the  medium.  Since  the  nature  of  the  Green’s 
function  changes  at  the  boundary  z  =  0  and  z  —  L,  the  energy  of  the  evanescent 
wave  would  return  back  into  the  media  pertinent  to  the  appropriate  region  of  the 
Green’s  function.  For  example,  at  z  <  0  no  energy  from  the  evanescing  waves  from 
reflection  would  go  into  the  medium,  nor  would  any  energy  from  the  evanescing  wave 
from  transmission  contribute  into  the  air. 

From  this  point  on,  we  will  ignore  the  effects  of  the  evanescent  waves  because 
we  will  assume  the  transmitter  and  receiver  are  sufficiently  far  away  from  this  region. 
But, we  can  see  that  in  a  half-space  *  <  0,  the  intensity  G*G  = 

becomes  smaller  as  z  gets  increasingly  negative.  We  expect  therefore  that  the  evanes¬ 
cent  contribution  will  be  small,  and  we  should  be  able  to  neglect  such  effects. 

F.  REDUCTION  OF  THE  INTEGRAL 

The  motivation  behind  the  coordinate  transformations  of  Equations  (IV.47), 
(IV. 48), (IV.51), (IV. 52)  is  twofold.  First,  they  will  allow  us  to  reduce  the  double 
integral  of  Equation  V.12  into  a  single  integral,  which  is  more  tractible.  But  we 
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also  do  this  transformation  in  conformity  to  Brekhovskikh’s  model.  This  approach 
will  allow  us  to  set  up  the  contour  integral  which  we  will  present  in  this  section  and 
perhaps  be  able  to  solve  in  the  future.  Since  the  method  for  which  this  transformation 
is  the  same  for  all  the  Green’s  functions,  we  will  demonstrate  it  for  the  case  of  the 
reflected  Green’s  function  in  this  section. 

The  symmetry  in  the  </>  direction  implies  a  switch  to  polar  coordinates: 

x  =  pcoscj)'  y  =  p  sin  <f>  (V.27) 


Substituting  this  equation  into  (V.12),  yields 


G(p,z,r0)  = 


i o  .to 


G(9,z,r o)ei(^sin0cos^co^+^sin0sin^'sin^  (  A  ]  cos  9  sin  6d9d(f> 

2  7 r 


Since 


27rJ0(fcpsin  9) 


We  are  left  with  a  Bessel  Transform: 


eikp  cos^-fl/)^ 


2 

2tt  J0(kp  sin  9)  cos  9  sin  9 d9  (V.28) 

Because  the  Hankel  functions  are  a  linear  combination  of  the  Bessel  series  of 
the  first  and  second  order,  Brekhovskikh  [Ref.  2]  writes  the  Bessel  Jo  in  the  following 
way: 


G(p,z,r0)  = 


G{9,z,r0) 


k 

2t r 


•Jo(kp  sin  9)  =  -  H^\kp  sin  9)  +  H^\kp  sin  9) 


r(2b 


(V.29) 


and,  using  the  Hankel  function  property 


H^\e-™x)  =  -H$\x), 


Equation  V.28  becomes 


k 


2  /*  7,  zoo 


l  r 


G(p,z,r0)  =  ^—  I  G(9,z,r0 W  (kp  sin  9)  +  H^’  (kp  sin  9) 
2tt  J0  2  1 


d2b 


cos  9  sin  9d9 
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k2 

Air 


zoo  p  Tj  —  ioo 

G(6,  z,  r0)  cos  9  sin  6H^\kp  sin  6,)d6*  +  /  G(0,  z,  r0)  cos  9  sin  6,h|)2') (kp  sin  #)d# 

J  o 


k2 
Ait 

The  reflection  Green’s  function  is  therefore 

G(p,2;,r0)  =  —  /  G(0, 2,r0)  cos6,sin0i/Q1)(A’psin^)d0 

4^  J-I+ioo 


G(0,  r0)  cos  9 sin  OH^^kp sin  6,)d6*  + 


(I- i°°) 


G(0,  z,  r0)  cos  6*  sin  HHq1'1  (kp sin  #)d0 


A’2 

47T 


m 


pikzo  cos  0  , - 

_ _ ikzyj  1— sin2  6 

2k  cos  9 


cos  9  sin  9H^  ( kp  sin  9)d9 


_  k2  r^°°  (1  -  n2)  +  (?r2  -  i)ei2fcLV/«2-sin2e 

^7r  J-^+ioo  (\/n2  —  sin2  9  +  cos 9)2  —  (a/ ri2  —  sin2  9  —  cos Qyiei2kL\/n2- sin2*? 

(V.30) 

pikzo  cos  6  , - 

X  — -  Q—ikzy/l—sin2  9  CQg  q  gin  qjjQ.)  ^  sin  q^q 

Zj  A;  COS  U 

G.  ANALYSIS 

Performing  the  contour  integration  of  Equation  V.30  requires  very  advanced 
mathematical  concepts.  One  has  to  approximate  the  integral  by  assuming  certain 
conditions  of  the  argument  of  the  integrand  —  such  as  wave  number  being  very  large. 
By  using  the  Method  of  Steepest  Decent  where  one  has  to  find  the  saddle  points  and 
determine  the  right  saddle  path,  an  analytical  approximate  solution  can  result  from 
Equation  V.30.  We  will  not  perform  this  integral  in  this  thesis,  but  we  will  analyze 
the  integral. 

First,  the  limits  of  integration  tell  us  that  the  summation  goes  from  all  the 
allowed  real  angles  to  the  complex  angles  (which  is  the  region  of  evanescent  waves). 
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Let  us  examine  the  latter  half  of  the  integral.  We  saw  earlier  that  Equation  V.3  is 
the  free  space  Green’s  function.  Analogously 


k 2 
An 


pikzQ  cos  6  I - 

— - -e~ikz V !— sin2  o  cos  6  sin  QH^ ( kp  sin  6)d6 

Zj  Hj  COS  \7 


(V.31) 


would  give  us  the  free  space  Green’s  function  (The  only  difference  being  that  it  is  left 
going  from  the  negative  sign  given  as  e~ikzV1-sin2  The  stunning  idea  is  that  the 
reflected  wave  is  spherical  in  nature  after  performing  the  Hankel  Transform.  This  is 
the  idea  behind  Huygen  wavelets  as  shown  in  Figure  11.  From  the  figure  you  can  see 
that  the  returning  waves  are  spherical,  and  they  superpose  with  each  other  to  create 
a  wave  front.  The  singularities  of  Equation  V.30  are  the  trapped  modes  of  either 
the  TE  or  TM  waves  in  the  media  (C.  Scandrett,  personal  communication,  December 
2006),  [Ref.  32],  This  approach  can  be  better  seen  when  converting  the  Cartesian 
form  into  the  polar  form  instead  of  the  spherical  form  (See  Appendix  G).  The  branch 
cuts  show  the  forbidden  regions,  and  these  are  the  regions  of  evanescing  waves  which 
Brekhovskikh  describes  in  meticulous  detail,  and  whch  require  a  contour  integration 
path  that  will  account  for  these  forbidden  regions  [Ref.  2]. 


©  W.  Fenclt  1998 

©  Prof.  T.  Mzouyhi  1998  Angle  of  Refraction:  14.5 0 


Figure  11.  Huygen  Wavelets.  Figure  was  taken  from  [Ref.  35] 


On  the  other  side  of  the  media  ( z  >  L),  we  can  also  use  Huygen  Wavelets. 
From  its  Cartesian  form,  we  can  demonstrate  see  this  effect  in  a  MATLAB  simula- 
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tion(see  Appendix  F).  In  Figure  12,  the  observer  is  on  the  plane  z  =  3m,  and  we  can 
see  these  wavelets. 


Transmission  v=500GHz,  Observation  point  at  z=3 


Figure  12.  Huygen  Wavelets.  Observation  at  z  —  3m 

It  maybe  beneficial  to  perform  the  integral  in  Equation  V.30  to  obtain  the 
Green’s  function  in  the  spatial  domain  because  it  shows  us  the  effects  of  all  these 
spherical  waves.  However,  this  thesis  will  focus  on  an  imaging  algorithm,  and  we  will 
work  in  the  spatial  frequency  domain  since  we  already  have  an  analytical  form.  We 
will  simulate  the  data  and  perform  the  transforms  numerically. 

H.  MATERIAL  ANALYSIS 

Listed  in  Table  II  are  some  common  barrier  materials.  We  will  examine  the 
reflection  from  the  wall.  As  expected,  the  type  of  material  will  greatly  affect  the 
amount  of  reflection  from  the  wall.  And  there  are  frequency  windows  at  which  we  get 
minimal  reflection.  Figures  13,  14,  and  15  show  the  reflectivity  of  materials  listed  in 
Table  II,  with  barrier  width  of  half  a  meter. 
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Frequency  vs  Reflection  from  the  dielectric 


Figure  13.  Reflection  (Cloth,  Balsa  Wood,  and  Wooden  Door) 

The  MATLAB  code  used  to  generate  these  results  uses  the  reflection  coefficient 
calculated  in  this  thesis.  The  code  is  modeled  for  a  plane  wave  parallel  to  the  surface, 
and  its  interaction  within  the  media. 

We  see  from  Figure  13  that  the  reflection  from  cloth  and  wood  are  very  small 
at  THz  frequencies  of  current  interest  and  so  the  material  becomes  almost  transparent 
to  these  waves. 

As  the  index  of  refraction  increases,  the  reflections  coming  back  from  the 
medium  become  more  significant  as  shown  in  Figure  14.  But  there  are  key  frequencies 
where  reflection  can  be  minimized.  For  the  case  of  Dry  Concrete,  a  common  material 
used  in  construction  of  most  walls,  the  material  essentially  becomes  transparent  at 
500  GHz.  But  a  potential  problem  relates  to  the  target  you  trying  to  image.  At 
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Frequency  vs  Reflection  from  the  dielectric 


Figure  14.  Reflection  (Plexiglas,  Paper,  Dry  Concrete) 

500  GHz,  plane  waves  will  penetrate  dry  concrete  with  little  loss  (at  one  meter  in 
depth).  If  you  want  to  get  information  on  people,  this  frequency  would  be  great.  If 
you  look  at  Figure  15,  muscle  is  not  very  reflective  at  this  frequency,  but  blood  is  as 
high  as  90%  reflective.  Thus,  gathering  data  scattered  from  people  at  this  frequency 
would  be  “ideal” .  In  general,  a  balance  must  be  achieved  between  wave  penetration 
of  the  dielectric  and  return  information  from  the  object — unless  the  object  is  metallic 
(in  which  case  the  only  concern  then  would  be  determining  the  best  penetrating 
frequency) . 
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Frequency  vs  Reflection  from  the  dielectric 


Figure  15.  Reflection  (Diamond,  Muscle,  Blood) 


Index  of  Refraction 

Cloth 

1.09 

Balsa  Wood 

1.14 

Wooden  Door 

1.41 

Plexiglas 

1.51 

Paper 

1.73  -  2.00 

Dry  Concrete 

2.12 

Diamond 

2.41 

Muscle  at  37  °C 

7.00 

Blood  at  37  °C 

7.62 

Table  II.  Typical  index  of  refraction  for  common  items 
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VI.  INTEGRAL  EQUATIONS 


A.  WAVE  EQUATION 

The  Electric  field  obeys  the  time  independent  Helmholtz  equation.  If  we  allow 
the  index  of  refraction  to  vary  so  slowly  with  position  that  it  is  essentially  constant 
over  distances  on  the  order  of  a  wavelength,  then  we  do  not  have  to  concern  ourselves 
with  the  coupling  of  Maxwell’s  Equations  [Ref.  25],  and  thus  we  may  examine  its 
scalar  form: 

V2£(r)  +  k2n2ohj(r)E(r)  =  0  (VI.l) 

We  can  re-write  this  equation  as 

V2E(r)  +  k2 E{r)  =  -  h2(nobj(r)2  -  1)  E(r)  (VI.2) 

' - -v - ' 

U 

where  U  denotes  the  “scattering  potential”  of  the  medium. 

V2£(r)  +  k2E(r )  =  -U(r)E(r)  (VI.3) 

The  the  Green’s  function  of  the  Helmholtz  equation  will  be  a  solution  of 

(V2  +  k2)G(r  -  7*0 )  =  -6(r  -  r0)  (VI.4) 

B.  LIPPMAN-SCHWINGER  EQUATION 

We  see  from  Figure  16  that  the  Green’s  function  is  divided  into  three  regions 
with  solutions  discussed  in  the  previous  chapter.  We  also  know  that  the  Green’s  func¬ 
tion  is  symmetric  so  that  the  response  at  r  due  to  a  concentrated  source  at  r0  is  the 
same  as  the  response  at  r0  due  to  a  concentrated  source  at  r.  We  will  use  these  ideas 
to  derive  the  Lippman-Schwinger’s  Equation  and  apply  the  first  Born  Approximation 
in  its  solution.  In  Equation  VI.3,  E(r)  is  the  total  field.  If  a  transmitter  is  in  region 
III  and  it  emanates  a  planar  wave  Einc(r )  =  Ae~lkz  into  the  media,  then  this  incident 
wave  will  be  planar  inside  the  medium  (region  II),  and  comes  out  as  a  planar  wave  on 
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the  other  side  into  region  I  as  E'inc(r)  =  tAe  lkz  where  t  is  the  transmittance.  This 
field  satisfies  the  homogeneous  Helmholtz  equation 

V2E'inc(r)  +  k2E’inc(r)  =  0  (VI.5) 

E(r)  can  be  expressed  as  the  sum  of  incident  field  and  the  scattered  field  so  that 

E(r)  =  E'inc(r)  +  E„  (r) 

So  where  Esc(r )  is  the  field  scattered  from  the  target.  Equation  VI. 3  can  be  therefore 
simplified  to 

V2Esc{r)  +  k2Esc(r)  =  -U(r)E(r)  (VI.6) 

Multiplying  Equation  VI. 4  by  Esc(r)  and  Equation  VI.6  by  G,  and  subtracting 
the  two  yields 

Esc(r)V2G(r  -  r0)  -  G(r  -  r0)V2EBC(r)  =  U{r)E(r)G(r  -  r0)  -  Eac(r)5(r  -  r0) 

(VI-7) 

Using  the  fact  that  the  Green’s  function  is  symmetric,  we  integrate  with  respect 
to  tq  throughout  the  entire  Volume  Vr  and  apply  Green’s  Second  Theorem  to  arrive 
at 


Esc(r)=  U(r0)E(r0)G(r0  -  r)d3r0 
JvR 


(VI.8) 


'Sr 


(■ E(r0)V°G(r0  -  r)  -  G(r0  -  r)V°Esc(r0))  ■  d S. 


R 


The  Sommerfeld  Radiation  condition  guarantees  us  that  the  fields  go  to  0  at  infinity, 
and  therefore,  as  the  surface  area  increases  with  the  radius  R,  the  area  integral 
vanishes.  Furthermore,  since  the  Scattering  potential  is  0  everywhere  outside  the 
scattering  object,  the  volume  integral  of  the  entire  volume  Vr  is  then  shrunk  down 
to  just  the  volume  of  the  scattering  object.  Therefore,  the  integral  reduces  to 


Esc(r)=  f  U(r0)E(r0)G(r0  -  r)d3r0  (VI. 9) 

Jv0 
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Figure  16.  Illustration  of  a  Scattering  Medium  occupying  Volume  Vr  bounded  by  a 
closed  Surface  Sr 

This  is  the  Lippman- Schwinger  Equation.  The  Green’s  function  is  different  for  each 
region,  depending  on  where  the  observation  point  is  located,  particularly  along  the 
z-axis: 


r 


G(r  -r0)  =  < 


Gi(r  -  r0) 
Gn(r  -  r0) 
Gin(r  -  r0) 


and  G/,  Gjj  and  G m  are  the  respective  Inverse 
(IV. 47),  (IV.3),  (IV. 51),  and  (IV.52). 


for  z  <  0 

for  0  <  z  <  L  (VI.  10) 

for  L  <  z 

Fourier  Transforms  of  Equations 
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C.  BORN  APPROXIMATION 


Substituting  Equation  VI. 9  into  Equation  Byields 


E  =  E'inc  +  /  U(r0)E(r0)G(r0  -  r) d3r0 
JVo 

If  we  formally  define  the  operator 

(AE)(r)  =  f  G(r o  —  r )E(r 0)d3r 0 
Jv0 

E  =  E'vru:  + 

Esc 

=  E',nc  +  A(E'^  +  AE) 

=  E'inc  +  AE'ine  +  A\E +  AE) 


(VI.11) 


Then  the  Born  approximation  of  Equation  VI.  11  becomes 

E  =  E'inc+[  U(r0)E'mc(r0)G(r0-r)d3r0  (VI.  12) 

J  Vo 

Esc=  [  U(r0)E'inc(r0)G(r0  -  r)d3r0  (VI.  13) 

JVo 


D.  THE  MODELING  SCHEME 

In  this  thesis,  we  will  only  account  for  single  scatter  events,  assuming  the  scat¬ 
tering  from  other  points  to  be  significantly  smaller  and  so  will  not  add  any  meaningful 
contributions.  We  shall  now  the  Scattering  Potential.  For  point  scatterers,  we  can 
easily  model  this  as  a  series  of  Dirac  Delta  Functions  so  that 


U(r)  =  Arf  (r  —  ri)  (VI.  14) 

i 

where  A  are  (complex- valued)  scatterer  amplitudes.  Substituting  Equation  VI.  14 
into  Equation  VI.  13  and  integrating,  we  get  a  summation  of  Green’s  functions.  Since 
the  receiver  will  be  in  region  III,  we  use  the  Transmitted  Green’s  function  in  region 


III. 


Esc{r)  -  J2AiAte~ikZiG^r  ~  r‘)  (VL15) 

i 
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VII. 


ALGORITHM 


A.  SPATIAL  DOMAIN  VS  FREQUENCY 

The  Green’s  function  forms  a  kernel  which  operates  on  the  Scattering  Potential 
and  the  Electric  held.  It  is  also  shift  invariant,  giving  us  the  special  characteristic  of 
the  convolution  integral.  Let  us  revisit  and  focus  on  the  volume  integral  of  Equation 


VI. 8 


Esc{r)  =  f  U(r0)E(r0)G(r0  -  r)d3r0  (VII.  1) 

JvR 


When  R  becomes  infinitely  large,  the  volume  integral  extends  over  all  space. 


OO 

U (:eo,  yo,  z0)E(x0 ,  yo,  z0)G(x-x0,y-y0,  £-z0)dxodyodz0  (VII. 2) 

—  OO 

It  is  the  scattering  potential  that  limits  this  integration  to  the  scattering  volume.  The 
first  Born  approximation  of  this  integral  yields 


Esc(x,  y,  z)  = 


U (x0,  2/o,  z0)E'inc(x0,  y0,  z0)G(x  -  x0,  y-y0,z-  z0)dx0dyodz0 

(VII. 3) 

Since  we  know  the  kernel  in  terms  of  its  spatial  frequencies,  we  perform  the 
convolution  integral  in  the  Fourier  domain.  We  reduce  the  3D  Fourier  Transform  to 
a  2D  Fourier  Transform  by  holding  z  constant.  The  idea  is  that  the  receiver  will 
be  placed  at  a  fixed  z  measurement.  The  2D  Fourier  Transform  of  this  convolution 
integral  will  result  in  the  2D  Fourier  Transform  of  the  Green’s  function,  and  a  2D 
Fourier  Transform  of  the  product  of  U  and  E'inc,  which  is  U  convolved  with  E[nc 
in  the  frequency  regime.  Since  our  scattering  model  is  a  summation  of  Dirac  Delta 
Functions,  the  2D  Fourier  Transform  of  Equation  VI.  15  with  respect  to  x  and  y  is 
simplified  to 


Esc(x,  y,  z)  = 
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E(fx,  fy)  =  Y,  AAte-^Gmif,,  /y)e-(2-/^+2-to)  (VII.5) 

i 

The  advantage  of  using  the  frequency  domain  is  that  we  deal  with  only  prod¬ 
ucts  in  MATLAB,  and  therefore,  no  fancy  coding  will  be  required  (simple  element-to- 
element  matrix  multiplications  are  enough).  We  will  also  let  the  computer  calculate 
the  2D  Inverse  Fourier  Transform  to  the  spatial  regime  and  compare  the  results. 


B.  MATLAB  CODE  IMPLEMENTATION 

Equation  VII.5  is  the  equation  we  shall  use  to  create  our  simulated  scattered 
field  data.  The  total  field  at  the  transmitter  will  be  composed  of  the  transmitted 
field  of  the  device,  the  reflected  fields  from  the  wall,  the  scattered  field,  and  the 
lateral  (evanescing)  field.  We  assume  that  the  evanescing  waves  will  be  minimal,  and 
that  we  are  purposely  staying  away  from  the  evanescing  regions.  Therefore,  the  only 
factors  we  account  for  are  the  transmitted,  reflected,  and  the  scattered  fields.  For 
our  simulation,  we  neglect  the  details  of  the  transmitted  fields  and  reflected  fields 
as  these  will  be  known  measurements  in  the  real  world.  For  image  processing,  we 
are  interested  in  the  scattered  field  information,  and  whether  Tikhonov  and/or  SVD 
methods  will  work  in  mitigating  noise  corruption. 

For  simulation  purpose,  we  reduce  the  propagation  number  so  that  the  effects 
of  the  point  scatterers  will  be  amplified.  The  limitation  of  this  model  is  that  it  is 
computation  intensive,  and  we  are  limited  to  the  sampling  frequency  and  bandwidth. 


64 


Point  Scatterers 


Figure  17.  Amplitude  of  field  measurements  due  to  5  Scattering  Points 

Figure  D  displays  the  noise- free  field  measurements,  collected  from  5  scattering 
points  measured  at  z  —  2m.  To  this  data,  we  add  Gaussian  noise,  and  we  will  apply 
the  method  of  Tikhonov  and  SVD  filtering  discussed  in  Chapter  2. 

The  code  is  presented  in  Appendix  D.  The  idea  is  to  do  all  the  work  in  the 
spatial  frequency  and  use  a  2-D  inverse  FFT  to  bring  the  results  back  to  the  spatial 
domain. 

C.  SIGNAL  TO  NOISE  RATIO 

The  additive  noise  is  Gaussian  noise  and  applied  in  the  frequency  domain. 
The  SNR  is  calculated  by  summing  up  the  intensities  of  each  element  of  the  matrix 
for  the  signal  and  the  noise,  and  taking  the  ratio  of  the  two.  This  provides  a  way  to 
gauge  our  data.  We  increased  the  Gaussian  noise  by  multiplying  it  by  a  noise  factor 
which  gives  us  the  desired  SNR  value. 
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Noisy  Data  with  SNR  of  -30.00 


Figure  18.  5  Scattering  Points  under  Gaussian  Noise  SNR=— 30dB 

D.  SPATIAL  FREQUENCY  DOMAIN 

Consider  the  case  of  a  two-point  scatterer.  For  one  scatterer,  the  frequency 
domain  follows  Figure  19.  The  second  point  scatterer,  looks  similar  to  the  first  one 
(see  Figure  20).  The  phase  shift  is  hard  to  observe. 

However,  when  we  add  the  two,  we  can  see  the  effects  of  phase  interference, 
(see  Figure  21,  and  upon  Fourier  inversion,  we  get  what  we  are  supposed  to  get  -  two 
scattering  points  (see  Figure  22). 

E.  TIKHONOV  REGULARIZATION 

From  Figure  18,  the  noise  has  almost  completely  masked  the  signal.  Now  we 
apply  Tikhonov  regularization,  and  the  we  will  choose  the  regularization  parameter  by 
trial  and  error.  We  know  that  as  the  regularization  parameter  decreases,  the  linear 
compact  operator  Ra  discussed  in  Chapter  2  will  head  towards  the  Least  Squares 
solution  to  the  problem  —  and  thereby  increase  the  condition  number.  Thus,  we  will 
see  that  the  recovered  image  will  become  more  unstable  (see  Figure  23). 
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Single  Point  Scatterer  in  Frequency  spectrum  (Transmission) 


Figure  19.  Single  Point  Scatterer  (Spatial  Frequency  Domain) 

This  figure  shows  the  ill-conditioned  nature  of  the  problem.  When  the  reg¬ 
ularization  parameter  is  increased  to  a  =  10,  we  see  how  well  the  noise  can  be 
mitigated  (see  Figure  24).  There  is  a  trade-off  between  a  large  regularization  param¬ 
eter  and  a  small  one.  A  larger  regularization  parameter  will  dominate  the  inverse 
operator  giving  you  optimization  problems,  while  a  small  one  will  give  you  a  least 
squares  solution  but  a  diverging  condition  number. 

F.  TRUNCATION  SVD 

We  now  apply  the  SVD-truncation  method.  At  first,  the  singular  values  were 
truncated  at  the  very  large  value  a  =  10000,  and  we  see  the  ill-conditioned  behavior 
(see  Figure  25). 

The  idea  is  to  compare  the  singular  values,  and  try  to  locate  where  the  rapid 
changes  of  singular  values  occur  and  truncate  at  that  point.  One  parameter  was 
determined  to  be  at  a  =  1,  and  the  image  is  pulled  out  of  the  noise  (see  Figure  26). 
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Figure  20.  Second  Point  Scatterer  (Spatial  Frequency  Domain) 

One  can  see  that  the  SVD  method  appears  to  do  a  visually  better  job  in 
retrieving  the  signal  from  the  noise,  but  this  may  not  be  generally  true. 
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Figure  21.  Two  Scattering  Points  (Spatial  Frequency  Domain) 
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Figure  22.  Two  Scattering  Points  (Spatial  Domain) 


Two  Point  Scatterer  in  spatial  domain  (Transmission) 
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Tikhonov  Regularization  alpha=0.003  and  SNR=-30.00 


Figure  23.  Tikhonov  Regularization  a  =  0.003  at  SNR=— 30  dB 


Tikhonov  Regularization  alpha=1 0.000  and  SNR=-30.00 


Figure  24.  Tikhonov  Regularization  a  =  10  at  SNR=— 30  dB 
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Figure  25.  TSVD  a  =  10000  at  SNR 


SVD  Truncation  w/  Singular  values  truncated  at  1 .00 
SNR=-30.00 


Figure  26.  TSVD  a  =  1  at  SNR=- 


SVD  T runcation  w/  Singular  values  truncated  at  1 0000.00 
SNR=-30.00 
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VIII. 


CONCLUSION 


Our  model  has  assumed  linearly  polarizable  media  with  unchanging  permit¬ 
tivity.  For  non-linear  materials,  the  polarization  vectors  will  have  to  be  accounted 
for  separately.  Though  tedious,  this  extension  would  provide  us  with  a  more  general 
model.  In  addition,  the  material  we  considered  is  not  magnetizable.  These  assump¬ 
tions  are  generally  good,  but  cannot  be  true  for  all  Through-the-Wall  cases.  It  is 
important  to  note  that  when  the  media  is  both  non-linear  and  magnetic,  the  Maxwell 
equations  will  be  coupled  and  the  scalar  theory  has  to  be  readjusted  to  fit  this  model. 
A  complementary  examination  of  common  wall  materials  will  be  required,  and  some 
of  the  future  work  is  to  see  how  well  the  model  works  in  “real-life”  applications,  and 
how  we  can  adjust  the  model  to  fit  specific  problems.  We  have  demonstrated  in  this 
thesis  that  the  scattered  field  can  be  modeled  as  simple  summation  of  Green’s  func¬ 
tions  because  we  applied  the  the  weak  scatterer  model.  When  the  Born  approximation 
is  not  valid,  we  may  have  to  consider  other  approximations,  and  we  may  also  want 
to  include  multiple  scattering  in  which  case  the  Lippman- Schwinger  equation  may  be 
extended  to  include  for  scatterers  from  nearby  points. 

The  Rytov  Approximation  assumes  that  the  wave  perturbation  is  caused  by 
objects  through  a  phase  change  in  the  incident  wave  along  the  propagation  path. 
It  is  derived  similarly  to  the  way  we  derived  the  Born  approximation,  except  the 
field  is  comprised  of  the  phases  as  functions  of  position.  It  is  said  that  the  Rytov 
approximation  is  a  more  accurate  approximation  than  the  standard  Born  model.  In 
future  work,  one  may  take  the  Lippman-Schwinger  Equation  and  apply  Rytov  and 
compare  to  the  Born  model.  It  may  be  worth  applying  this  to  “real”  world  problems 
along  with  the  Born  model  for  future  work. 

We  discussed  non-conductive  dielectrics.  For  conductive  media,  the  attenua¬ 
tion  maybe  too  great  for  examining  the  scattered  field.  The  attenuation  of  a  signal 
through  the  media  and  back  maybe  masked  by  the  Gaussian  noise  and  the  data  may 
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be  too  noisy  for  any  reasonable  image  reconstruction  algorithm  to  work.  For  direct 
application  to  real  world  scenarios,  the  conductive  case  may  be  unlikely  to  be  en¬ 
countered.  For  imaging  behind  walls,  most  construction  materials  are  not  completely 
made  of  conductive  materials.  However,  steel  bars  may  pose  a  problem  with  back 
scatter  inside  the  media,  which  we  would  have  to  deal  with  separately. 

We  also  ignored  dispersive  effects  that  need  to  be  examined  for  various  mate¬ 
rials  for  which  a  priori  knowledge  of  the  material  characteristics  must  be  known  in 
order  to  make  the  appropriate  approximations  for  the  theory.  It  will  be  of  interest  in 
future  work  to  make  appropriate  approximations  and  account  for  dispersion  by  using 
the  Method  of  Stationary  Phase.  Such  approximate  analytical  solutions  can  provide 
some  understanding  about  how  the  waves  behave  in  transit  through  the  medium,  but 
also  can  be  much  more  efficient  and  accurate  in  the  image  reconstruction  stage. 

The  Tikhonov  regularization  parameter  was  found  by  trial  and  error.  There 
are  much  more  effective  ways  to  determine  the  regularization  parameter  such  as  the 
L-curve  method  [Ref.  30].  The  Bayesian  statistical  approach  [Ref.  30]  to  the  problem 
can  give  us  better  image  reconstruction  and  might  be  fruitful  for  future  research. 

Our  code  implementation  has  some  limitations.  Since  MATLAB  can  be  very 
memory  hungry,  the  bandwidth  examined  could  not  be  very  wide,  nor  can  one  increase 
the  sampling  frequencies  very  high  without  sacrificing  time  (and  patience).  A  more 
efficient  computation  scheme  will  be  required.  In  this  thesis,  I  have  used  short  band- 
widths  with  sampling  frequencies  at  Nyquist  frequencies  to  avoid  aliasing.  A  dual 
core  2Ghz  with  6GB  of  RAM  would  take  5  —  15  minutes  to  do  one  full  computations 
for  a  reasonable  frequency  bandwidth  of  20m-1  for  fx  and  fy.  This  poses  a  problem  in 
computation  and  analysis  on  systems  requiring  higher  bandwidths.  Efficient  coding 
and  higher  computing  power  will  be  required  for  analysis  of  larger  bandwidths. 
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APPENDIX  A.  THE  RIEMANN  LEBESGUE 

LEMMA 


lim  /  K(X,  s)C  sm(ns)ds  =  0 


(A.l) 


Through  integration  by  parts,  we  get  the  following  and  therefore  conclude  that  as 
71  — »  CO,  the  integral  will  then  goes  to  0. 


lim  — 

n — »oo 


K(X,  s)  cos(ns)  |b  1 


71 


a  n 


9K(X,s)  ,  w  n 

.  - - - cos  (7is)ds  =  0 

2  '  a.s 


(A.2) 
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APPENDIX  B.  CURRENT  IMAGING 
METHODS  (ANL) 


At  Argonne  National  Laboratory,  A  receiver  and  transmitter  via  waveguides 
are  used  for  imaging.  The  wave  guides  are  shown  in  the  following  diagram.  The 
configuration  is  set  to  send  electromagnetic  waves  from  10  Ghz  to  200  Ghz,  but  it 
also  has  the  ability  to  go  to  higher  frequencies  if  necessary.  Both  the  transmitter  and 
receiver  are  of  the  same  type  and  design,  and  they  are  adjacent  to  each  other  while 
in  radar  mode  (i.e.  data  received  from  scattered  fields).  In  the  radar  mode,  both  the 
transmitter  and  receiver  are  moved  simultaneously  in  space  in  sequence. 


Figure  27.  Waveguide  for  transmitter.  Separate  waveguide  for  receiver  not  presented 
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APPENDIX  C.  SINGULAR  VALUE 
DECOMPOSITION 


Let  A  be  a  linear  operator  that  maps  X  — »  Y,  and  to  be  more  specific,  let 
us  say  that  X  C  ]Rn  and  Y  C  ]Rm.  So  that  A  will  have  dimension  m  x  n.  In  this 
section  we  will  discuss  ATA  ( n  x  n )  and  AAT  (m  x  m ). 

The  eigen  equation  for  ATA  is 

ATAx.i  =  A  iXi  (C.l) 

where  A i  is  the  eigenvalue  corresponding  to  the  vector  Xi.  Because  AT  A  is  symmetric, 
the  eigenvectors  are  orthogonal  [Ref.  29].  Similarly,  the  eigen  equation  for  AAT  is 

AATyj  =  ipjVj  (C.2) 

where  ifjj  is  the  eigenvalue  corresponding  to  y ;.  AAT  is  also  a  symmetric  matrix  and 
therefore  the  eigenvectors  yi are  orthogonal. 

Consider  the  eigenvector  Xi  of  AT A,  with  nonzero  eigenvalue  Aj.  Then  Axi  ^ 
0  because  Ax.j  is  an  eigenvector  of  AAT: 

(AAt)Ax2  =  A(ATA)Xi  =  A  i(Axi)  (C.3) 

The  same  is  true  for  an  eigenvector  y  of  AAT  with  nonzero  eigenvalue  since  ATy  is 
an  eigenvector  of  ATA. 

It  follows  that  the  non-zero  eigenvalues  of  ATA  must  also  be  eigenvalues  of 
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AAt  and  vice  versa.  If  the  rank  of  ATA  is  r,  then  this  means  that 

Ai  =  i/’i 

A2  =  '02 

Ar  =  0r 

Ar+i  0  'i/v+i  0 


(C.4) 

(C.5) 

(C.6) 


Then  equation  C.5  becomes 

ATyk  =  akxk  Axk  =  akyk  (C.7) 

Since  the  eigenvalues  of  \k  =  0  and  t/jk  —  0  for  k  >  r,  this  means  that 


Axk  =  0  for  k  =  r  +  1, . . .  n 
ATyk  =  0  for  k  =  r  +  1, . . .  m 

If  we  take  Axk  =  okyk ,  and  multiply  equation  C.7  by  xT  on  both  sides,  and 
use  equation  C.8,  we  see  that 

r 

A  =  aiyiXl  (C-8) 

i=  1 
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and  similarly, 


r 

AT  =  ^2  (J>x>yl  (C-9) 

i= 1 

The  orthonormal  vectors  x  are  called  right  singluar  vectors  and  y  are  called 
left  singular  vector.  In  matrix  form,  we  can  write  equation  C.8  as 


or 

A  =  YSXt  (C.ll) 

where  S'  is  a  diagonal  matrix  of  the  singular  values,  and  Y  and  X  are  composed  of 
the  eigenvectors  corresponding  to  those  singular  values. 

The  inverse  is  easily  calculated  by  the  following: 

A~l  =  ( YSX T)"1 

Since  XT  and  Y  are  orthogonal  matrices  due  to  the  fact  that  the  eigenvectors  are 
orthogonal,  we  see  the  following: 

A~l  =  (XT)-1S~1Y-1 

A-1  =  X5_1yT  (C.12) 

Combined  with  the  Least  Square  Method,  this  result  can  be  used  to  get  us  a  solution 
to  the  system  of  equation.  As  you  can  see  the  inverse  mapping  provides  us  with  ^ 
for  its  diagonal  matrix,  and  this  becomes  a  problem  as  the  singular  values  get  smaller 
and  smaller.  Since  the  non-zero  eigenvalues  of  a  Hilbert  Space  operator  are  isolated 
and  form  a  sequence  that  converges  to  0  [Ref.  28],  the  norm  of  inverse  mapping 
A-1  whose  singular  values  of  -,  will  increase  indefinitely.  Therefore,  regularization 
methods  need  to  be  introduced  to  cure  this  effect. 
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APPENDIX  D.  MATLAB  CODE  FOR 
SCATTERING  POINTS 


clear ; 


°/0Parameters 
nu=500e9 ; 
lambda=3 . 0e8/ nu ; 
7ok=2*pi/lambda; 

k=3000 ; 
n=2 . 4 ; 

L=0 . 5 ; 
z=2 ; 

A=l; 

signal=0 ; 
noise=0 ; 


*/.  Frequency  of  the  wave  propagation. 

%  In  terms  of  wave  length. 

*/.  Magnitude  of  the  propagation  vector. 

%  Material  index  of  refraction. 

*/.  Thickness  of  the  wall. 

I  Measurement  plane  at  z=  constant.  Adjust  accordingly. 
*/.  depending  on  position  of  the  receiver. 

I  Amplitude  of  the  incident  plane  wave. 


7„  Defining  Coordinate  Systems  in  the  spatial  frequency  domain. 
ii=20;  7o  spatial  frequencies  for  meshgrid  input. 

fb=ii;  7o  frequency  band. 

7of  n=f  b ; 


fn=2*fb;  7o  Nyquist  rate. 

j=l/fn;  7o  Separation  must  be  at  least  this  distance  to  prevent  aliasing, 

[f x , fy]  =meshgrid(-ii  :  j  :  ii  , -ii  :  j  :  ii)  ;  °/0  Creating  coordinate  system. 


7o  Defining  the  Green's  function. 

x0=0 ; y0=0 ; z0=-3 ;  7o  Set  this  position  accordingly  to  the  problem. 

sN=sqrt (k~2- (2*pi*f x) . “2- (2*pi*f y) . ~2) ;  7o  Look  in  chapter  3  for  reference. 

sNm=sqrt ( (n*k) ~2- (2*pi*fx) . ~2- (2*pi*f y) . ~2) ;  7o  Look  in  chapter  3  for  reference. 
num=4*sN . *sNm*exp (-i*sN*L) ;  °/0  Numerator  of  the  equation. 

7.  Denominator  of  the  equation. 

den=(sNm+sN) . ~2 . *exp(-i*sNm*L) - (sNm-sN) . ~2 . *exp(i*sNm*L) ; 

X0=(exp(-i*2*pi*fx*x0) . *exp(-i*2*pi*fy*y0) . *exp(-i*sN*zO) . *exp(i*sN*z) ) ,/(2*sN)  ; 


G=(num.  /den)  .  *X0;  °/0  Transmission  Green's  function. 
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%%r/.r/.r/.mmmrm  Synthetic  Data  Generation 


0/0  /o  /o  /o/o 


%  transmittance  for  a  plane  wave  through  a  given  media 
t=4*n*exp(-i*k*L) . /((n+1) ~2*exp (-i*n*k*L) - (n-1) ~2*exp (i*n*k*L) ) ; 


7.  Scattering  Points  in  the  spatial  frequency  domain. 

xl=-l;yl=-l;zl=-l; 

x2=-l ; y2=l ; z2=-3 ; 

x3=0 ; y3=0 ; z3=-5 ; 

x4=l;  y4=l ; z4=-l . 3 ; 

x5=l;  y5=-l ; z5=-2 . 2 ; 

7„  Incident  field  after  the  penetration  of  the  media 
El=t*A*exp(i*k*zl) ; 

E2=t*A*exp(i*k*z2) ; 

E3=t*A*exp(i*k*z3) ; 

E4=t*A*exp(i*k*z4) ; 

E5=t*A*exp(i*k*z5) ; 


7.  Positions  of  the  scatter  points  in  frequency  domain  for  the  given 
7.  measurement  plane  at  z=constant 

7o  Since  the  scatter  points  are  in  terms  of  dirac  delta  functions, 

7.  the  field  measurement  would  become  the  summation  of  the  Green’s  function 
7.  and  the  incident  wave  after  performing  the  Born  Approximation 


Xl=(exp(-i*2*pi*fy*xl) . *exp(-i*2*pi*fx*yl) . *exp(-i*sN*zl) . *exp(i*sN*z) ) ,/(2*sN) ; 
X2=(exp(-i*2*pi*fy*x2) . *exp(-i*2*pi*fx*y2) . *exp(-i*sN*z2) . *exp(i*sN*z) ) ,/(2*sN) ; 
X3=(exp(-i*2*pi*fy*x3) . *exp(-i*2*pi*fx*y3) . *exp(-i*sN*z3) . *exp(i*sN*z) ) ,/(2*sN) ; 
X4=(exp(-i*2*pi*fy*x4) . *exp(-i*2*pi*fx*y4) . *exp(-i*sN*z4) . *exp(i*sN*z) ) ,/(2*sN) ; 
X5=(exp(-i*2*pi*fy*x5) . *exp(-i*2*pi*fx*y5) . *exp(-i*sN*z5) . *exp(i*sN*z) ) ,/(2*sN) ; 

7.  Generating  Data 

h=(num . / den) . * (E1*X1+E2*X2+E3*X3+E4*X4+E5*X5) ; 


7oSetting  matrix  dimension 
[M,N] =size(h) ; 

h=h(l  :M,  1  :N)  ;  7oChange  matrix  dimension  if  warranted. 

7oConverting  bin  numbers  to  real  axis  numbers  in  spatial  domain  (meters)  . 
70Range  has  to  be  up  to  Nyquist  frequency. 
ax=linspace (-fn/2 , fn/2 , length(h) ) ; 
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n=(randn(M,N)+i*randn(M,N) )  ;  Ttioise  in  the  frequency  domain. 


°/0  Set  the  signal-to-noise  (in  dB) 
snr=-35 ; 


°/0Signal  to  Noise  ratio 

SIG=abs (h . *conj (h) ) ; 

N0I=abs (n . *conj (n) ) ; 
signal=sum(sum(SIG) ) ; 
noise=sum(sum(NOI) ) ; 
f actor=sqrt ( (signal/noise) *10" 
n=n*f actor ; 

d=h+n ; 


-snr/10) ) ; 

°/0data  in  the  frequency  domain. 


nnnilllinnniniT/1  End  of  (noisy)  Data  Generation70%%%%r/o0/o0/o0/o7.7o7.7.7.7.7.7.7o 


figure(l)  7o  data  corrupted  with  noise  in  spatial  domain. 

dl=ifft2(d);  %  data  in  the  spatial  domain. 

dl=if ftshift (dl) ; 

imagesc(ax,ax,abs(dl)) ; 

colormap  (1-gray); 

axis ([-2  2-22]) 

title  (sprintf  (’  Noisy  Data  with  SNR  of  702 . 2f 1  ,  snr) ) 

7.  Tikhonov  Regularization 

alpha=3;  %  Regularization  parameter  for  Tikhonov. 

R=inv(alpha*eye (size (G ’ *G) ) +G’ *G) *G’ ; 

figure (2) 

p=R.*d; 

p=ifft2(p);  7«  data  in  the  spatial  domain. 

p=if ftshift (p) ; 

imagesc(ax,ax,abs(p)) ; 

axis ([-2  2-22]) 

colormap  (1-gray); 
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title  (sprintf  ( 1  Tikhonov  Regularization  alpha=°/02 . 3f  and  SNR=°/02 . 2f  alpha,  snr)  ) 


°/0  SVD  Truncation  Method 

alpha=.l;  %  Parameter  for  Truncation  of  singular  values, 

[u, s , v] =svd(G) ; 

7„  inverse  of  the  singluar  value  matrix. 
sinv=zeros (N,M) ; 
for  l=l:rank(s) 

if  1/s (1 , 1) <alpha 

sinv(l , l)=l/s(l , 1) ; 

else 

sinv(l,l)=0; 

end 

end 


e=(v*sinv*u’ ) .  *d; 

e=ifft2(e);  %  data  in  the  spatial  domain. 

e=if ftshift (e) ; 


figure (3) 

imagesc(ax,ax,abs(e)) ; 
colormap  (1-gray); 
axis ([-2  2-22]) 

title  (sprintf  ( 1 TSVD  w/  singular  values  cutoff  at  %3.2f  \n  SNR=°/02 . 2f  1  ,  alpha,  snr) ) 
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APPENDIX  E.  MATLAB  CODE  FOR 
MATERIAL  REFLECTION 


clear ; 

"^Calculating  the  Reflection  for  plane  waves  in  various  surface 
0/0Equations  used  from  Chapter  3 

°/0Creating  various  index  of  refraction  for  different  materials 

70Cloth,  Balsa  Wood,  Wooden  Door,  Plexiglas,  Paper,  Dry  Concrete,  Diamond, 

°/„Muscle ,  Blood 

nm= [1 . 09  1.14  1.41  1.51  1.73  2.12  2.41  7.00  7.62]; 

°/0Thickness  of  the  material 
L=0 . 5 ; 

c=3.0e8;  °/0speed  of  light 


“/oCloth,  Balsa  Wood,  Wooden  Door 
figure (1) 
for  j = 1 : 3 

nu=linspace (20e9 , lel2 , 200) ; 
k=2*pi*nu . /c ; 

num=(  (l+nm(  j ) )  *  (l-nm(  j ) )  +  (nm(  j  )  +  l)  *  (nm(  j )  -1)  *exp  (i*2*k*nm(  j )  *L) )  ; 
den=( (nm( j )+l) * (nm( j )+l) +(nm( j ) -1) * (l-nm( j ) ) *exp (i*2*k*nm( j ) *L) ) ; 
r=num . / den ; 

R=r . *conj (r) ; 
plot (nu,R) ; 

text (nu(30* ( j  —  1 ) +8) ,R(30* ( j — 1 ) +8) , sprintf  ( ’n=  %2 . 2f 1 ,nm( j ) ) ) ; 
hold  on 


end 

hold  off 

xlabel (’ frequency  (\nu)  Hz’); 
ylabel( ’Reflection’ ) ; 

title( ’Frequency  vs  Reflection  from  the  dielectric’) 
text  (lell ,.  115 ,  sprintf  (’ L  =  °/01.2f  m’,L)); 

°/„axis (  [0 .  lell  lOell  0  1]); 
clear  figure 


°/0Plexiglass  ,  Paper ,  Dry  Concrete 
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figure (2) 
for  j=l:3 

nu=linspace (20e9 , lel2 , 100)  ; 
k=2*pi*nu . / c ; 

num=( (l+nm( j+3) ) * (l-nm( j+3) )+ (nm( j+3)+l) * (nm( j+3) -1) *exp(i*2*k*nm( j+3) *L) ) ; 
den=( (nm( j+3)+l) * (nm( j+3)+l)+ (nm( j+3) -1) * (l-nm( j+3) ) *exp(i*2*k*nm( j+3) *L) ) ; 
r=num . / den ; 

R=r . *conj (r) ; 
plot (nu,R) ! 

text  (nu(30*  (j-l)  +  10)  ,R(30*  ( j — 1) +10)  ,  sprintf  ( ’n=  °/„2 . 2f  ’  ,nm(j+3) )  )  ; 
hold  on 


end 

hold  off 

xlabel (’ frequency  (\nu)  Hz’); 
ylabel( ’Reflection’ ) ; 

text  (lell 4,  sprintf  (’ L  =  °/01.2f  m’,L)); 

°/„axis (  [0 .  lell  lOell  0  1]); 

title( ’Frequency  vs  Reflection  from  the  dielectric’) 


°/0Thickness  of  the  material 
L=0 . 01 


°/0Diamond,  Muscle,  Blood 
figure (3) 
for  j = 1 : 3 

nu=linspace (20e9 , lel2 , 100)  ; 
k=2*pi*nu . / c ; 

num=( (l+nm( j+6) ) * (l-nm( j+6) )  +  (nm( j  +6) +1 ) * (nm( j+6) -1) *exp(i*2*k*nm( j+6) *L) ) ; 
den=( (nm( j  +6 ) + 1 ) * (nm( j+6)  +  1)  +  (nm( j+6) -1) * (l-nm( j+6) ) *exp(i*2*k*nm( j+6) *L) ) ; 
r=num . / den ; 

R=r . *conj (r) ; 
plot (nu,R) ; 

text  (nu(30*  (j-l)  +  10)  ,R(30*  ( j — 1) +10)  ,  sprintf  ( ’n=  °/„2 . 2f  ’  ,nm( j+6) ) )  ; 

°/„text  (nu(6*  ( j+2) )  ,R(  (j+5) )  ,  sprintf  (  ’n=  °/„l .  2f  ’  ,nm(  j+6)  )  )  ; 
hold  on 


end 


hold  off 

xlabel (’ frequency  (\nu)  Hz’); 
ylabel( ’Reflection’ ) ; 

title( ’Frequency  vs  Reflection  from  the  dielectric’) 
text  (lell ,  .  95 ,  sprintf  (  ’  L  =  °/01.2f  m’,L)); 
axis(  [0 . lell  lOell  0  1]); 
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APPENDIX  F.  MATLAB  CODE  FOR 
TRANSMITTED  WAVE  FRONT 


clear ; 

yoThis  code  graphs  the  transmission  of  the  waves 
“/Observation  of  the  fields  at  (0,0,3) 

/(Parameters 

nu=500e9;  %  Frequency  of  the  wave  propagation. 

lambda=3 . 0e8/nu;  %  In  terms  of  wave  length. 

k=2*pi/lambda;  %  Magnitude  of  the  propagation  vector. 

°/„k=3000 ; 

n=2.4;  %  Material  index  of  refraction. 

L=0.5;  %  Thickness  of  the  wall. 

z=3;  l  Measurement  plane  at  z=  constant.  Adjust  accordingly. 

I  depending  on  position  of  the  receiver. 

signal=0 ; 
noise=0 ; 

%  Defining  Coordinate  Systems  in  the  spatial  frequency  domain. 
ii=1000;  %  spatial  frequencies  for  meshgrid  input. 

j=-89; 

°/0j=l/fn;  %  Separation  must  be  at  least  this  distance  to  prevent  aliasing. 

[f x , fy] =meshgrid(-ii : j : ii , -ii : j : ii) ;  %  Creating  coordinate  system. 

[x,y] =meshgrid(-ii :j:ii,-ii:j:ii); 

%  Defining  the  Green's  function. 

x0=0 ; y0=0 ; z0=-3 ;  %  Set  this  position  accordingly  to  the  problem. 

sN=sqrt (k~2- (2*pi*f x)  .  "2- (2*pi*f y)  .  ~2)  ;  °/0  Look  in  chapter  3  for  reference. 

sNm=sqrt ( (n*k) ~2- (2*pi*fx) . ~2- (2*pi*f y) . ~2) ;  %  Look  in  chapter  3  for  reference. 
num=4*sN . *sNm*exp (-i*sN*L) ;  %  Numerator  of  the  equation. 

%  Denominator  of  the  equation. 

den=(sNm+sN) . ~2 . *exp(-i*sNm*L) - (sNm-sN) . "2 . *exp(i*sNm*L) ; 

X0=(exp(-i*2*pi*fx*x0) . *exp(-i*2*pi*fy*y0) . *exp(-i*sN*zO) . *exp(i*sN*z) ) ,/(2*sN) ; 

G=(num. /den) . *X0;  %  Transmission  Green's  function. 

figure (1) 
imagesc (abs (G) ) 
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colormap (1-gray) 

title (’  Transmission  \nu=500GHz,  Observation  point  at  z=3') 
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APPENDIX  G.  GREEN’S  FUNCTION  POLAR 
COORDINATE  TRANSFORMATION 


If  we  make  the  following  coordinate  transformation: 


fl  +  fy=P 2 


fx  =  P  cos  (f)  fv  =  p  sin  < 


2,2  2 

x  +  y  =  r 


x  =  r  cos  0  y  =  r  sin  0 

Then  the  reflection  Green’s  function  of  the  Equation  IV. 47  can  be  written  as 


G(r,  9,  z,  r*0)  =  2ir  /  pdp  ~ ' — -x 


x 


J0(27rpr)  k2(n2  -  l)(ei2LV/^Z^V  _  i)e-i(^)Vk2-4^ 

2  y/k_47r2p2 

1 


(a/ n2k 2  —  47rp2  +  \Jk2  —  Air2  pi2)2  —  (y/n2k2  —  Air2p2  —  ^Jk2  —  Air2  p2)2el2L^n2k2  47r2p2 
Let 


2irp  nku  nkdu 

u  =  — r-  P=  dp  = 


2vr 


2vr 


G  = 1 1  £ )  *v  - 11 1 


uduJo(nkru)  ^el2nkL%/1  u2  —  1^  e  l(^z  zo)kVT 


fc\/l  — 


2/,.2 


n*u 


x 


fc2(?r\/l  —  u2  +  \/l  —  n2u2)2  —  k2{ns/ 1  —  u2  —  \/l  —  n2n2)2ez2Lnkv/l 


G  = 


x 


kn2(n2  —  1)  [°°  MduJo(nkru)  ^el2nkLx/l  u2  —  1^  e  zo)kVi  n2u2 


47T 


\/l  — 


2n,2 


n*u 


(ny/ 1  —  u2  +  Vl  —  n2u2)2  —  (ns/ 1  —  u2  —  \/T—  n2u2)2el2Lnk^~ 
let  7l  =  n-\/l  —  u2  and  B  =  \/l  —  n2u2  0  =  nkL\Jl  —  u2 
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Where  are  the  zeroes  of  the  denominator  of  the  integrand? 
where  does 


(. A  +  B )2 

A2  +  B2  +  2  AB 
2  AB 
A2  +  B 2 


Denominator  vanishes  when 


(A  -  Bfe 

(. A 2  +  B2  - 
ei2e  -  1 

ei2Q  _|_  l 

1  2  sin  0 

2  cos  0 


i20 

-  2AB)ei2e 

g^®  _  g 

e*e  +  e*e 
i  tan  0 


i  t&n(nkLy/l  —  u 2)  = 


_  2ny  1— u2Vl— n2z 


n2  +  l— 2n2u2 


Look  at  RHS  Numerator  is  real  and  positive  for  0  <  u  <  1  jn 
it  then  becomes  imaginary  and  positive  for  1/n  <  u  <  1 
it  then  becomes  real  and  negative  for  u  >  1 

Denominator  goes  to  zero  at  u  —  u*  =  yj 


Is 


So  RHS, 


u*  <  1  for  n  >  1 


?r2  +  1 
2n2 
n2  +  1 


2?r2 


<  1 


<  1 


n 2  >  1 


u*  >  1/n? 


n2  +  1  ?  1 
2?r  n: 


n2  +  1  >  2 


n 


>  1 


1  *  1 

-  <  u  <  1 

n 

N 

D 


0  <  u  <  1/n 
1/n  <  u  <  u* 
u*  <  u  <  1 
u  >  1 


real, positive 
positive 

imaginary, positive 
positive 
imag,pos 
neg 

real,neg 

neg 


real  >  0 
imag  >  0 
imag  <  0 
real  >  0 
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Look  at  the  LHS  ita,n(nkLy/l  —  u 2) 


0  <  u  <  1  LHS  imaginary,  either  positive  or  negative 


tan(nr)  =  itanh(aj) 

u  >  1  i  ■  ita.nh(nkLy/u2  —  1)  =  —  tanh (nkLy/ u2  —  1)  always  negative,  real 

The  conclusion  only  possible  to  have  a  root  for  1/n  <u  <  1 
What  about  the  end  points? 

Always  a  root  at  u  =  1  0  =  0 

For  u  —  1/n 


i  t&nLikLx  1 - -)  =  itan(kLVn2  —  1) 

V  nz 


2ny/l  —  \/n2iy/n2u2  —  1  2  iyjn2u2  —  1 

n2  -  1  >/n2  -  1 

as  u  — »  1/n  from  above 

tan(fcL\/?r2  —  1)  =  0 

possible  to  have  a  root  if  kLy/n 2  —  1  =  imr  Always  a  root  at  n  =  1. 

Root  at  u  =  1/n  if  kLy/n2  —  1  =  mir  where  m  is  an  integer. 

Range  of  possible  roots  1/n  <  u  <  1 
How  many  roots  are  there? 

The  “real”  equation  we  need  to  examine  is 

tan(nfcLVl  -  n2)  =  - y— - — ^ - 

nA  +  1  —  zn^ir 

To  analyze  this  equation,  lets  change  variables  let 

y  =  nkLy/ 1  —  u 2 


u2  =  1  — 


nkL 


Limits:  as  u  —>  1,  y  — ►  0.  u  — >  1/n,  y  — >  nkL^J  1  —  ^  =  kLy/n 2  —  1 
Gall  ymax  =  kLy/n-  1  ,  0  ^  y  ^  j/max 


„2  +  1  -  2n2(1  _  ( J_)2) 


/„2  _  1  _  y2 

kL  Y  1  (kL)2 


— (?r2  —  1)  + 


2y2 

(kL)2 


2 y  y/ (kL)2(n2  -  1)  -  y2 


(kL)2  —  (n2  -  1)  + 


2y2 

(kL)2 
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RHS 


2y\/ymax-y2 

ymax2+2y2 


Thus,  we  see  roots  of 


tan  y 


2y  V  i/max-y2 

-'/max  +  2’/2 


(Note  compilation  with  the  aid  of  Professor  James  Luscombe) 


This  transcendental  equation  is  analogous  to  the  same  transcendental  equation 
which  solves  the  TE  modes  for  the  Dielectric  Slabs  [Ref.  32] .  The  intersection  of  the 
two  curves  between  the  left  and  right  hand  side  of  the  equation  are  the  values  for 
which  propagation  in  the  slab  is  possible.  The  branch  cuts  of  the  integral  represents 
the  forbidden  regions  where  propagation  is  not  possible  and  where  the  waves  begin 
to  evanesce.  Brekhovskikh  and  Kong  work  out  in  detail  of  these  effects  [Ref.  2]  [Ref. 
32],  The  problem  with  this  integral  is  the  nature  of  the  Bessel  function.  As  we 
integrate  to  infinity,  the  Bessel  function  also  head  towards  infinity,  and  therefore, 
there  is  no  way  to  close  the  contour  integral.  More  clever  methods  are  required,  and 
Brekhovskikh  model  provides  a  way  to  perform  the  contour  and  to  implement  the 
contour.  Nevertheless,  we  can  see  that  the  poles  are  the  trapped  modes  in  the 
media. 
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